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ABSTRACT 


Jr 

A  system  of  computer  programs  has  been  developed  to  predict  flows 
around  and  In  the  near  wake  behind  an  axlsymmetric  body  including  a  pro¬ 
peller  at  the  tail  of  the  vehicle.  The  system  is  capable  of  calculating 
laminar  and/or  turbulent,  attached  and  separated  flows.  Interactive 
boundary-layer  theory  is  used  to  predict  the  flow  on  the  nonseparating 
region  over  the  body  and  the  full  Navier-Stokes  equations  are  solved  for 
the  separated,  propeller  and  near  wake  flow  regions. 

Details  of  the  theory  and  results  are  included  in  a  separate  en¬ 
gineering  report.  The  programs,  subroutines,  variables  in  common  and 
input  and  output  data  are  described.  Listings  of  the  program  code.  In¬ 
put  and  output  data  for  a  sample  case  are  Included. 


FORWARD 


The  following  system  of  computer  programs  has  been  developed  for 
use  on  the  IBM  370,  model  158  computer,  for  CMS  under  VM.  For  this 
reason  the  exec  procedures  required  to  execute  this  system  have  also 
been  included  to  aid  in  the  conversion  of  the  system  in  its  present  form 
to  a  form  compatible  with  that  in  use  at  a  receiving  facility. 
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INTRODUCTION 


The  purpose  of  this  system  of  programs  is  the  prediction  of 
viscous  flow  over  and  in  the  near  wake  behind  axisymmetric  bodies, 
particularly  submarine-like  bodies.  The  effects  of  an  operating  pro¬ 
peller  and  separated  flow  have  been  included  in  the  procedure.  The 
system  procedure  is  basically  the  execution  of  four  segments.  The  first 
is  the  development  of  a  transformed  or  nonseparating  body  which  allows  a 
boundary- layer  solution  to  be  calculated.  The  second  segment  is  the 
solution  of  the  interactive  boundary-layer  over  the  body.  The  third 
segment  is  the  development  of  the  perturbed  flowfield,  used  in  the 
fourth  and  final  segment,  the  solution  of  the  afterbody  and  near  wake 
with  the  propeller. 
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DESCRIPTION  OF  SYSTEM 


The  system  is  composed  of  a  number  of  programs  which  serve  not 

i 

only  to  obtain  a  solution,  but  to  as  great  a  degree  as  possible  automate 
the  calculation  of  the  solution.  The  first  segment  of  the  system  is 
composed  of  support  programs  included  specifically  to  perform  this  task. 
The  general  body  coordinates  are  input  and  these  programs  calculate  a 
transformed  set  of  coordinates  that  allow  a  nonseparating  boundary- layer 
solution  to  be  obtained.  The  basis  of  the  transformation  is  in  the 
nature  of  the  flow  over  the  afterbody.  The  original  coordinates  are 
tested  to  determine  if  the  flow  separates  over  the  afterbody.  The  test 
is  based  on  the  inviscid  pressure  distribution  over  the  body  and  was 
empirically  developed  from  repeated  trials  on  a  number  of  test  configura¬ 
tions.  From  these  trials  it  was  determined  that  the  boundary  layer  on 
the  afterbody  would  not  remain  attached  if  the  inviscid  pressure  coef¬ 
ficient  exceeded  a  value  of  0.2,  Therefore,  an  inviscid  calculation  is 
performed  to  give  the  pressures  over  the  afterbody  and,  if  the  criterion 
for  separation  is  met,  the  body  is  modified  until  a  nonseparating  con¬ 
figuration  is  found. 

The  next  segment  of  the  system  is  an  interactive  boundary-layer 
solution  about  the  transformed  body.  The  interactive  solution  employs 
the  VPI&SU  Incompressible  Boundary-Layer  Program  and  the  Hess  Potential 
Flow  Program  to  obtain  a  viscous/inviscid  iterated  solution.  The  itera¬ 
tion  procedure  is  as  follows:  A  boundary-layer  solution  is  first  ob¬ 
tained  by  matching  the  boundary-layer  edge  pressure  with  a  "first  guess" 
inviscid  pressure  distribution.  Coordinates  of  the  displacement  body 
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that  result  from  the  boundary-! ayer  solution  are  placed  in  a  separate 
coordinate  data  set.  The  inviscid  pressure  distribution  is  then  found 
about  the  displaced  body  by  the  potential  flow  program  and  becomes  the 
next  guess  for  the  pressure  distribution.  This  guess  is  then  used  by 
the  boundary-layer  program  to  yield  another,  corresponding  displaced 
body  and  the  sequence  is  continued.  The  iterative  process  that  is  de¬ 
veloped  is  continued  for  several  iterations,  the  number  of  which  is  de¬ 
pendent  on  the  shape  of  the  afterbody  and  the  accuracy  of  the  initial 
guess.  With  the  current  system,  the  first  guess  is  the  inviscid  pres¬ 
sure  distribution  about  the  original  body  obtained  from  the  potential 
flow  code.  It  has  been  found  that  the  use  of  this  initial  guess  requires 
that  only  one  iteration  was  needed  when  solving  for  the  flows  about 
bodies  with  straight  tails  and  from  4-6  iterations  were  needed  for 
straight-tapered  or  curved  tails. 

The  elliptic  nature  of  the  Navier-Stokes  equations  requires  that 

full  (surrounding),  boundary  conditions  be  supplied,  as  well  as  a  first 

guess  of  the  flowfield  in  the  region  treated.  The  downstream  boundary 

lies  in  the  wake  of  the  body  and  the  corresponding  boundary  conditions 
2  2  2  2 

are  that  a  c/3*  and  3  i|i/3x  go  to  zero.  Since  the  region  that  is 
treated  in  the  Navier-Stokes  program,  in  this  case  a  "box"  of  grid 
points,  encompasses  a  much  greater  area  than  that  of  the  viscous  region 
alone,  the  full  viscous-inviscid  velocity  profile  is  required  at  the  up¬ 
stream  edge.  The  interactive  boundary- layer  solution  provides  conditions 
for  the  viscous  region  of  the  upstream  boundary,  and  conditions  in  the 
inviscid  region  are  obtained  by  employing  the  potential  flow  program  to 
calculate  the  outer  flowfield  about  the  final  displaced  body.  At  the 
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joining  of  the  two  velocity  profiles,  the  viscous  and  inviscid  regions, 
a  composite  expansion  is  used  to  ensure  a  continuous  transition.  The 
initial  flowfield  is  thus  an  inviscid  velocity  field  and  is  used  inside 
the  Navier-Stokes  program  to  generate  the  initial  stream  function 
throughout  the  Navier-Stokes  region.  The  upper  boundary  conditions  are 
simply  the  inviscid  velocities  obtained  from  the  perturbed  flowfield. 

The  final  segment  of  the  system  is  the  solution  of  the  propeller, 
separated  and  near  wake  regions  of  the  flow  which  is  accomplished 
through  the  use  of  the  APLNS  program.  The  Navier-Stokes  program  iterates 
on  the  stream  function,  which,  on  the  first  iteration,  is  the  stream 
function  about  the  inviscid,  perturbed  body  described  above.  An 
Alternating  Difference  Implicit  method  (ADI)  is  used  for  calculating 
new  stream  functions  in  the  Navier-Stokes  (N-S)  region.  After  NTS  time 
steps,  the  solution  is  assumed  to  reach  convergence  and  the  results  are 
printed.  For  many  high  speed  problems  the  time  step  (At)  must  be 
particularly  small,  and  hence  the  number  of  time  steps  required  for 
convergence  may  be  considerable. 

There  is  no  final  convergence  criterion  other  than  the  specified 
total  number  of  time  steps  (NTS). 

The  logical  flow  of  the  system  may  be  followed  in  Fig.  1  . 

Figure  2  shows  the  corresponding  sequence  of  program  usage  in  the 
system.  This  user's  guide  contains  a  brief  overview  of  the  system 
augmented  by  a  more  complete  and  detailed  description  of  the  system  and 
its  component  programs.  The  system  is  tied  together  by  WAKE  EXEC,  an 
exec  procedure  which  executes  the  segment  procedures.  See  Appendix  A, 
Table  1.1  for  a  listing  of  WAKE  EXEC. 
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Figure  I.  LOGICAL  FLOW  OF  SYSTEM 


Figure  2.  MECHANICAL  FLOW  OF  SYSTEM 


DESCRIPTION  OF  PROCEDURE  CDCP 

CDCP  EXEC  is  an  exec  procedure  which  ties  together  the  operation 
of  the  support  programs:  DATADEV,  DATAMA,  and  HESS. 

Programs  DATADEV  and  DATAMA,  controlled  by  CDCP  EXEC,  generate 
a  complete  set  of  data  files  with  which  the  iterative  procedure  can  be 
executed.  The  procedure  includes  the  development  and  the  checking  of 
body  geometries  until  a  nonseparating  flow  is  found.  The  development 
of  the  nonseparating  flow  involves:  (1)  the  incremental  extension  of 
a  trailing  edge  fairing  into  the  wake  (performed  in  DATADEV),  (2)  the 
testing  of  that  fairing  in  HESS,  and  (3)  the  check  of  that  fairing  in 
DATAMA.  CDCP  EXEC  first  executes  DATADEV,  then  HESSBLI  and  then  DATAMA 
execs.  When  DATAMA  has  completed  executing,  the  data  file  SEPARA  DATA 
is  checked  and  the  flag  ICP  is  evaluated.  IF  ICP  is  0,  a  successful 
afterbody  has  been  found;  if  1,  the  iterative  procedure  is  continued. 

If  a  nonseparating  fairing  has  been  found,  data  sets  are  written  and 
control  shifts  to  BLINT  EXEC,  but  if  not,  the  fairing  is  extended  further 
into  the  wake,  and  the  result  is  retested.  This  sequence  is  continued 
until  a  successful  afterbody  is  generated  or  an  iteration  limit  is 
reached.  Because  of  imperfections  in  the  fitting  of  the  faired  after¬ 
body  to  the  original  body,  as  stated  above,  it  is  suggested  that  the 
user  carefully  examine  the  coordinate  sets  and  Cp  distribution  that  are 
output  by  CDCP  and  smooth  where  needed,  especially  in  the  trailing  edge 
region  and  near  the  body/fairing  joint. 

Use  of  CDCP  EXEC  (ARGUMENT) 

(ARGUMENT)  is  the  maximun  number  of  iterations  that  the  sequence 
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will  complete  to  find  a  satisfactory  fairing.  This  is  put  in  as  a  safety 
feature  to  prevent  execution  from  continuing  indefinitely  due  to  some 
unforeseen  body  characteristic.  A  listing  of  CDCP  EXEC  can  be  found  in 
Appendix  A,  Table  2.1. 

Program  DATADEV 

Description  of  Program  DATADEV.  DATADEV  EXEC  is  the  EXEC  program 
which  executes  DATADEV  FORTRAN.  DATADEV  FORTRAN  generates  data  sets  for 
the  HESS  potential  flow  code.  On  the  first  pass  DATADEV  simply  writes 
out  a  full  HESS  data  set  using  the  original  body  coordinates.  On  suc¬ 
cessive  passes,  DATADEV  generates  the  afterbody  which  extends  into  the 
wake  and  further  up  the  body.  The  extension  up  the  body  has  been  found 
to  hasten  afterbody  development.  The  fairing  is  concave  upward,  begin¬ 
ning  at  the  point  where  Cp  changes  sign  from  negative  to  positive  (see 
Fig.  3  ).  After  developing  the  afterbody,  new  HESS  data  sets  are  written 
out,  and  a  new  data  set  with  iteration  information  is  written.  A  flow 
chart  of  DATADEV  is  found  in  Fig.  4  . 

See  Appendix  A,  Tables  2.2  and  2.3  for  listings  of  DATADEV  EXEC 
and  DATAMA  FORTRAN. 

Description  of  Variables  in  Common. 


BLOCK 

VARIABLE 

DESCRIPTION  OF  USE 

CORFIT 

C 

Slope  of  body  at  the 
beginning  of  the 
Navier-Stokes  region. 

CORFIT 

DX 

Constant  stepsize 
used  in  re-positioning 
afterbody  coordinates. 

REPOSITION 
END  COORDINATES 


Figure  4.  LOGICAL  FLOW  OF  DATADEV 
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COORD 

KTE 

Station  number  of 
trailing  edge  on 
original  body. 

COORD 

KTE1 

Station  number  of 
trailing  edge  of 
faired  (extended) 
body. 

CORF IT 

NNN 

Station  number  of 
beginning  of  re¬ 
positioned  coordinates. 

COORD 

NOPRFS 

Number  of  profiles 
in  the  part  of  the 
Navier-Stokes  region 
on  the  body. 

COORD 

X 

Original  x  coordinates 
on  first  pass,  re¬ 
positioned  coordinates 
on  subsequent  passes. 

COORD 

Y 

Original  y  coordinates 
on  first  pass,  re¬ 
positioned  coordinates 
on  subsequent  passes. 

CORFU 

XI 

Axial  distance  from 
the  beginning  of  the 
fairing  to  KTE 1 . 

CORFU 

Y1 

Radial  distance  from 
the  beginning  of  the 
fairing  to  the  axis. 

COORD 

XI 

Array  of  x  coordinates 
with  afterbody  fairing. 

COORD 

YI 

Array  of  y  coordinates 
with  afterbody  fairing. 

Description  of  Subroutines. 

MAIN.  MAIN  reads  in  the  body  coordinates.  On  the  first  pass, 
the  original  coordinates  are  read.  On  successive  passes  the  original 
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coordinates  with  the  coordinates  of  the  afterbody  repositioned  to  give  a 
constant  x  stepsize  are  read.  MAIN  then  calls  LINDEV  and  HESSPA 
to  generate  and  write  out  new  coordinates  with  afterbody. 

The  variables  listed  below  are  used  in  MAIN.  A  description  can 
be  found  in  VARIABLES  IN  COMMON:  KWAKE,  KTE ,  KTE1 ,  X,  Y. 

Subroutine  LINDEV.  LINDEV  develops  the  trial  afterbody 
fairings  and  reads  in  the  separation  information,  including  the  itera¬ 
tion  number.  The  routine  then  repositions  the  coordinates  near  the 
trailing  edge  to  give  a  constant  stepsize  and  on  subsequent  iterations, 
extends  the  fairing  into  the  wake  LWAKE  points  and  steps  back  the  be¬ 
ginning  of  the  fairing  IWAKE  points.  This  set  of  x-coordinates  is  used 
in  developing  the  fairing  which  is  done  in  FIT. 

Description  of  Variables. 

DX  X  stepsize  in  afterbody  for  repositioned  coordinates.  Set 

at  one  percent  of  body  length. 

KK  Number  of  repositioned  coordinates. 

LWAKE  Number  of  points  in  wake  used  in  extension  of  afterbody, 

LWAKE  =  IWAKE*5 

MM  Number  of  body  stations  with  repositioned  coordinates. 

Subroutine  FIT.  Subroutine  FIT  generates  a  curve  to  fit  the 
boundary  conditions  given  by  LINDEV  in  the  form  of  XI,  X2,  Yl,  and  C. 

C  is  the  body  slope  at  XI,  Yl;  the  slope  at  the  trailing  edge  is  assumed 
to  be  zero.  The  fit  is  based  on  the  equations. 

Y^  =  CaY(K)  +  aY(K2)  +  aY(K3)  +  ...  +  aY(KN)], 

N  =  1,  KK,  KK  =  No.  of  points  on  curve  (1) 


(2) 
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K  =  ly  C  (3) 

DELY  is  initially  set  to  0.0  and  successively  increased  by  DY. 

K  is  calculated  according  to  Eq.  3.  If  Y  from  Eqs.  1  and  2  is  less 
than  Y 1 ,  DELY  is  increased  and  the  calculation  is  redone  until  the  con¬ 
vergence  criterion  (YK-Y1 )  =  0.0001  is  attained.  The  Y's  based  on  the 
last  K  are  used  in  the  afterbody  coordinate  sets. 

Description  of  Variables. 

DELX  X  stepsize  (=  dx) 

DELY  Y  stepsize 

K  Constant  multiplier  in  numerical  curve  fit 

L  Station  number  of  the  beginning  of  the  afterbody 

Y  Working  array  of  Y  coordinates  of  afterbody 

Other  variables  used  in  FIT  and  described  in  Variables  in  Common 

are  C,  Yl,  XI,  YI,  and  DX. 

Subroutine  HESSPA.  HESSPA  writes  out  a  complete  HESS 
data  set  using  the  body  coordinate  arrays  XI  and  YI  which  include  the 
afterbody.  The  routine  first  reads  in  the  HESS  parameters  and  includes 
the  changed  parameters,  if  any,  in  the  written  set. 

Description  of  Input  and  Output. 

Input  Data.  Three  data  sets  are  input  into  DATADEV  for 
execution.  They  are:  COORD  DATA,  SEPARA  DATA,  and  HESS  PARMDATA. 

COORD  DATA  contains  the  original  body  coordinates  together  with 
parameters  and  flags. 

Card  1  KWAKE  (13) 
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Card  2 


Cards 

3-KTE+3 


KWAKE 


KTE,  NOPRFS, 
KTE1  (lx, 313) 


NOPRFS 


Flag  for  including  wake  points  in  solu¬ 
tion:  0  no  wake  points 
1  use  wake  points 

(must  have  data  set  of  wake  points  if 
KWAKE  =  1) 


Station  number  of  trailing  edge  on  original 
body  (first  pass  original  T.  E.,  suc¬ 
cessive  passes  repositioned  T.  E.) 

Number  of  profiles  on  the  body  in  the 
APLNS  box.  Maximum  of  15  is  the  recom¬ 
mended  value 

Station  number  of  trailing  edge  on  faired 


(2F12.6) 


X  x  coordinate  of  the  original  body  geometry 

Y  y  coordinate  of  the  original  body  geometry 

See  Appendix  B,  Table  1.1  for  a  sample  COORD  DATA  file. 

SEPARA  DATA  contains  information  on  the  status  of  the  iterations 


and  flags  on  separation. 


Card  1 


Card  2 


I  WAKE 


JSTA,  KTE 


Iteration  counter,  on  first  iteration  has 
a  value  of  0 
(lx, 313) 

Station  number  of  start  of  faired  body. 
For  initial  data  set  can  be  0 


i 
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' 

KTE 

Described  above 

Card 

3 

ICP 

(13) 

I  CP 

Separation  flag,  0  for  no  separation; 

1  for  separation.  Set  to  0  initially. 

See  Appendix  B,  Table  1 

.2  for  a  sample  SEPARA  DATA  data  set. 

HESS 

PARMDATA  contains 

the  HESS  execution  parameters.  A  further 

description 

can  be  found  in  the  HESS  user’s  guide. 

Card 

1 

HEDR,  CASE 

(15A4,2x,2A4) 

HEDR 

Alphanumeric  array  with  the  case  title 

CASE 

Alphanumeric  array  with  the  case  number 

Card 

2 

ICARD2 

(13) 

ICARD2 

HESS  execution  flags,  for  body  only  solu 

tion  set  at  1 ,0,1 

Card 

3 

ICARD3 

(15) 

ICARD3 

HESS  execution  flags,  set  to  0  (see  HESS 

user's  guide) 

Card 

4 

ICARD4,  NN 

(lx, 14, 15) 

ICARD4 

HESS  execution  flags,  for  body  only  solu 

tion  set  to  1 ,0,0,1 

NN 

Number  of  body  stations.  Set  to  KTE 

initially  and  will  be  changed  to  KTE1 

for  separating  case 

Card 

5 

IBDN 

(9x.Il) 

IBDN 

Body  number,  for  body  only  solution  set 

to  1 

See  Appendix  B,  Table  1.3  for  a  sample  HESS  PARMDATA  data  set. 
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Output  Data.  The  output  files  created  by  DATADEV  are: 


SEPARA  DATA,  COORD  DATA,  and  HESSBLI  DATA. 

SEPARA  DATA  Described  above.  Rewritten  for  use  by 


COORD  DATA 


HESSBLI  DATA 


Card  1-5 


DATAMA 

Described  above.  Values  of  KTE1  and 
original  coordinates  changed  to  reflect 
fairing  extension  (KTE1)  and  repositioning 
of  original  body  coordinates 
Complete  data  set  for  HESS  execution, 
including  fairing  when  one  is  used 
(second  and  subsequent  iterations) 
Described  above  in  HESS  PARMDATA 


Card  6  -  KTE1  XI,  YI 
XI  ,YI 


(2F12.6) 

Body  coordinates 


See  Appendix  B,  Table  1.4  for  a  sample  HESSBLI  DATA  data  set. 


Program  Datama 

Description  of  Program  Datama.  DATAMA  EXEC  executes  program 
DATAMA.  DATAMA  evaluates  the  pressure  distribution  developed  by  HESS 
and  decides  if  the  geometry  will  produce  separating  flow.  The  program 
also  interpolates  the  Cp  data  set  at  on-coordinate  points  for  use  in  the 
HESS/ICBLINT  iterative  scheme.  If  the  developed  afterbody  success¬ 
fully  prevents  separation,  the  full  HESS/ICBLINT  data  sets  are 
written.  An  option  is  included  to  add  wake  points  for  a  wake  solution 
from  HESS/ICBLINT.  A  flow  chart  of  DATAMA  is  shown  in  Fig.  5  . 

See  Appendix  A,  Tables  2.4  and  2.5  for  listings  of  DATAMA  EXEC  and 


DATAMA  FORTRAN. 
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Description 

Variable 

DYDX 

JSTA 


I  CP 

UMAX 


JWAKE 

KTE 

KTE1 

KWAKE 

X 

XI 

XLOCAT 

Y 

YI 


of  Variables  in  Common. 

Description  of  Variable 
Slope  of  straight  taper  approximation  of 
tail  section  (without  afterbody)  for 
use  in  APLNS  geometry  data 
Station  near  trailing  edge  at  which  the 
Cp  changes  sign  (-  to  +),  which  determines 
the  beginning  of  the  fairing  extension. 
Separation  flag;  1  if  flow  is  separating, 

0  if  otherwise 

Maximum  number  of  coordinate  stations. 
Values  in  coordinate  array  including  wake 
points 

Number  of  wake  points  at  which  a  solu¬ 
tion  is  to  be  found 
Station  number  of  original  body 
Station  number  of  trailing  edge  of  body 
with  fairing 

Wake  flag;  1  if  wake  points  are  to  be 

used,  0  if  otherwise 

Array  of  faired  body  coordinates 

Array  of  original  body  coordinates 

X  coord  of  beginning  of  fairing  X(JSTA) 

Array  of  final  y  body  coordinates  with 

fairing 

Array  of  original  body  coordinates 


14 


YWAKE  Array  of  Y  values  of  wake  coordinates. 

Description  of  Subroutines. 

MAIN.  MAIN  reads  in  wake  coordinates,  if  used,  and  calls 
subroutines  CPDATA,  CURVE,  and  READAT.  The  last  two  routines  are  called 
only  if  ICR  =  0. 

Subroutine  CPDATA.  Subroutine  CPDATA  is  the  working 
section  of  the  program.  Pressure  coefficients  are  read  in  from  the 
previous  HESS  execution  and,  since  these  coefficients  are  given  at  off- 
coordinate  points,  the  data  are  interpolated  to  on-body  coordinate 
points.  These  Cp  data  are  examined  to  find  locations  of  sign  changes, 
negative  minimums  and  magnitudes  over  the  fairing.  If  the  fairing  Cp 
data  go  higher  than  0.2,  the  body  is  assumed  to  produce  separating  flow. 
If  wake  points  are  to  be  used,  they  are  added  to  the  coordinate  array 
and  given  a  constant  Cp  value  of  the  last  negative  minimum. 

Description  of  Variables. 

CP  Array  of  interpolated  Cp  values 

F  Array  of  Cp  values  for  off-body  points 

ICPAR  Station  number  of  negative  minimum  Cp 

used  in  giving  initial  values  to  wake 

C  data 
P 

I WAKE  Iteration  counter 

LWAKE  Number  of  points  in  wake  used  in  fairing 

extension  =  IWAKE*5 

Subroutine  CURVE.  CURVE  is  called  when  ICP  =  0.  CURVE 


reads  in  the  original  body  coordinates  and  finds  the  slope  of  the  body 
at  XI(KSTA).  The  procedure  used  defines  the  straight  tapered  tail  for 
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the  APLNS  geometry.  CURVE  also  prints  the  necessary  data  for  DEVELOP 
to  calculate  the  inviscid,  perturbed  flowfield. 

Subroutine  READAT.  Subroutine  READAT  reads  in  the  HESS/ 
ICBLINT  execution  parameters  and  writes  out  the  complete  working  data 
sets.  For  coordinates  the  routine  uses  the  last  (now  nonseparating) 
faired  set  developed  in  DATADEV  found  in  HESSBLI  DATA.  The  routine 
also  writes  out  the  stations  at  which  complete  solutions  in  ICBL  are  to 
be  obtained  (IPR  array)  and  the  stations  at  which  solutions  are  to  be 
written  out  (IPRINT  array).  For  simplicity  these  arrays  are  equal,  the 
stations  used  are  dependent  on  UMAX. 

Description  of  Variables. 

K  Station  number  at  which  IPRINT  array 

begins  to  be  incremented  by  1. 

Description  of  Input  and  Output. 

Input  Data.  There  are  seven  data  sets  used  by  DATAMA:  WAKE 
COORDATA,  ICBLINT  PARMDATA,  HESSBLI  DATA,  SEPARA  DATA,  SUMDS  DATA, 

HESSBLI  OUTPUT  and  COORD  DATA.  Of  these  data  sets,  HESSBLI  DATA, 
SEPARA  DATA,  and  COORD  DATA  have  been  described  above  in  INPUT  and  OUTPUT 
for  DATADEV. 

WAKE  COORDATA.  WAKE  COORDATA  contains  the  coordinates  of 
wake  points  if  a  wake  solution  is  desired.  The  input  data  are  as  follows: 
Card  1  JWAKE  (13) 

JWAKE  JWAKE  is  the  number  of  wake  points  in  the 

arrays  XWAKE  and  YWAKE 

Card  2-KWAKE+l  XWAKE  (F12.6)  KWAKE  is  the  number  of  points 

in  wake  coordinate  data  set 
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XWAKE  The  array  of  values  of  coordinates  in  the 

array  of  wake  points. 

See  Appendix  B,  Table  1.5  for  a  sample  WAKE  COORDATA  data  set. 
ICBLINT  PARMDATA.  ICBLINT  PARMDATA  is  the  data  set  of 
HESS/ICBL  execution  parameters.  The  data  set  contains  all  the  necessary 
ICBLINT/HESS  data  except  for  the  axial,  radial,  and  surface  coordinates. 
The  data  are  updated  in  DATAMA  to  include  any  changes  made  in  the  coordi¬ 
nates.  For  a  complete  description  of  ICBLINT  PARMDATA  see  the  ICBLINT 
user's  guide,  "Description  of  Input,"  under  SUB  DATA.  SUB  DATA  contains 
in  Cards  1-20  the  set  ICBLINT  PARMDATA. 

SUMDS  DATA.  SUMDS  DATA  is  the  data  set  output  by  HESS 
giving  the  surface  distance  to  each  of  the  coordinates  which  becomes  the 
XSTA  array  used  in  ICBLINT.  The  data  set  is  read  in  by  READAT  and  written 
on  SUB  DATA. 

Card  1-NN  XSTA  (12x,Fl2.6) 

See  Appendix  B,  Table  1.6  for  a  sample  SUMDS  DATA  data  set. 

HESSBLI  OUTPUT.  HESSBLI  OUTPUT  is  the  off -coordinate  Cp 
data  from  HESS. 

Cards  1-NN  XW  +  XTO  «-  F(3F12.6) 

XW  XW  is  the  off-coordinate  point  at  which 

the  C  data  have  been  calculated  in  HESS 
P 

XTO  XTO  is  the  on-coordinate  point  to  which 

the  Cp  data  will  be  interpolated 
F  F  is  the  off -coordinate  Cp  value 

Output  Data.  DATAMA  outputs  five  data  files:  SUB  DATA, 
HESSCP  OUTPUT,  SEPARA  DATA,  DATAMAKE  DATA,  and  INTER  DATA. 
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SUB  DATA.  SUB  DATA  is  the  complete  data  file  for  HESS/ICBL 
execution  with  all  parameters  and  coordinate  points  (see  INPUT  FOR 
ICBLINT) . 

HESSCP  OUTPUT.  HESSCP  OUTPUT  is  the  final,  interpolated 
data  set  of  Cp  values  for  use  in  HESS/ICBL  along  with  SUB  DATA.  The 
data  contain  the  on-coordinate  Cp  values.  HESSCP  OUTPUT  is  described 
in  the  ICBLINT  user’s  guide,  "Description  of  Input." 

See  Appendix  B,  Table  1.7  for  a  sample  HESSCP  OUTPUT  data  set. 

SEPARA  DATA.  SEPARA  DATA,  described  earlier,  is  changed  in 
DATAMA,  and  output  again  with  the  new  information. 

DATAMAKE  DATA.  DATAMAKE  DATA  is  the  data  set  for  execution 
of  BLDMKER  FORTRAN  and  is  simply  the  number  of  points  in  the  final  body 
coordinate  set. 

Card  1  NN  (13) 

NN  NN  is  the  number  of  body  points 

See  Appendix  B,  Table  1.8  for  a  sample  DATAMAKE  DATA  data  set. 

INTER  DATA.  INTER  DATA  is  the  data  set  for  DATFIELD  to  use 
in  developing  the  inviscid  flowfield. 

Card  1  IHESSI  (511.15) 

IHESSI  Counter  for  the  number  of  repetitions 

completed  through  DEVELOP 
Card  2  XLOCAT,DX,DY  (3F12.6) 

XLOCAT  Location  of  the  upstream  edge  of  Navier- 

Stokes  region 

DX  x  stepsize  for  placement  of  off-body 

points 


DESCRIPTION  OF  PROCEDURE  BUNT 


BLINT  EXEC  is  the  controlling  exec  procedure  for  the  HESS/ICBLINT 
system  of  codes.  The  Iterative  scheme  involves  the  execution  of  these 
programs,  together  with  the  support  program  BLDMKER  in  a  loop.  The  num¬ 
ber  of  loops  is  user  determined,  since  an  accurate  criterion  for  con¬ 
vergence  is  difficult  to  define.  In  using  the  method,  it  has  been  found 
that  generally  convergence  Is  obtained  after  the  sixth  iteration.  For 
simple  geometries,  fewer  than  six  iterations  are  required. 

The  HESS/ICBLINT  iterative  sequence  is  the  prime  data  developer 
of  the  APLNS  program.  In  the  present  form  the  sequence  uses  the  HESS 
potential  flow  program  in  very  nearly  its  original  form,  the  exceptions 
being  the  output  of  off-body  point  solutions,  pressure  distributions,  and 
surface  coordinates.  Program  ICBLINT  has  been  modified  to  both  accept 
and  employ  In  calculations,  Cp  values  from  the  inviscid  solution,  calcu¬ 
lation  of  wake  properties  and  a  restart  capability.  All  of  these  are 
used  in  the  viscous/inviscld,  iterated  solution.  The  concepts  and  theory 
of  this  iterative  procedure  have  been  described,  and  the  following  is  a 
description  of  the  code  mechanics. 

Figures  6  and  7  show  the  operation  of  BLINT.  The  procedure 
begins  with  the  execution  of  ICBLINT  to  yield  the  initial  displaced  body 
B00YC00R  DATA,  which  Is  run  by  HESS.  HESS  output  is  the  Cp  distribution 
on  the  body.  The  Cp  data  set  (HESSBLI  OUTPUT)  is  manipulated  by  BLDMKER 
to  yield  HESSCP  OUTPUT  which,  together  with  the  original  coordinates  and 
the  displaced  body  coordinates,  forms  the  Input  for  the  second  iteration. 
In  the  second  and  subsequent  Iterations,  the  original  body  and  the  dis¬ 
placed  body  coordinates  are  averaged  together  (weighting  factor  of  0.5) 
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OBTAIN  BOUNDARY 
LAYER  SOLUTION 


Figure  6.  LOGICAL  FLOW  OF  HESS/ICBLINT 


Figure  7.  MECHANICAL  FLOW  OF  HESS  M  CBL I  NT 
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in  ICBLINT.  Each  run  yields  a  displaced  body  (again  B0DYC00R  DATA)  which 
is  used  as  HESS  input.  The  procedure  then  repeats  itself  for  the  number 
of  iterations  desired. 

The  procedure  for  the  wake  calculations  enables  the  user  to  obtain 
iterated  wake  solutions  for  nonseparating,  no-propeller  cases.  The  ex¬ 
tension  into  the  wake  may  be  as  far  as  desired  within  the  array  dimensions 
which  will  be  discussed  later.  Data  set  WAKE  C00RDATA  containing  the  wake 
points  has  been  described  earlier. 

The  restart  capability  is  a  method  of  saving  execution  time  on  the 
second  and  subsequent  executions  of  ICBLINT.  There  is  virtually  no  inter¬ 
action  on  the  forward  region  of  most  bodies  and  therefore  no  reason  to  redo 
boundary-layer  calculations  in  this  region.  Therefore,  a  tape  is  written 
at  the  restart  station,  and  in  subsequent  iterations  the  tape  is  read 
and  the  solution  begun  at  the  restart  station.  The  actual  location  of 
this  station  is  dependent  on  geometry  and  must  be  determined  initially  by 
the  user. 

See  Appendix  A,  Table  3.1  for  a  listing  of  BUNT  EXEC. 

Program  HESS 

The  HESS  potential  flow  code  is  controlled  by  HESSBLI  EXEC.  Only 
three  modifications  have  been  made  to  the  original  version  of  the  Douglas 
computer  program  and  are  as  follows: 

1.  The  pressure  coefficients  are  written  out  on  a  separate 
tape  (HESSBLI  OUTPUT). 

2.  The  complete  inviscid,  off -body  flowfield  is  written  on  a 
separate  tape,  and  added  on  during  each  execution.  Details 
of  this  modification  are  found  in  the  description  of  program 
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DEVELOP,  where  the  nodif ication  is  used. 

3.  The  surface  distance  from  the  leading  edge  (SUMDS  DATA)  is 
written  out  on  separate  tape  for  use  in  ICBLINT  input  data 
sets. 

The  complete  HESS  user's  guide  is  sent  under  separate  cover, 
describing  the  code  and  is  accurate  except  for  the  above  modifications. 

A  listing  of  HESSBLI  EXEC  is  given  in  Table  3.2  of  Appendix  A. 

Program  BLDMKER 

Description  of  Program  BLDMKER.  BLDMKER  is  an  interpolating  pro¬ 
gram  which  performs  the  same  function  as  Program  DATAMA  in  routine  CPDATA. 
BLDMKER  interpolates  the  Cp  data  output  by  HESS  in  HESSBLI  OUTPUT  and 
finds  the  Cp  at  on-coordinate  positions,  creating  HESSCP  OUTPUT. 

See  Appendix  A,  Table  3.3  and  Table  3.4  for  listings  of  BLDMKER 
EXEC  and  BLDMKER  FORTRAN. 

Description  of  Input  for  BLDMKER.  BLDMKER  requires  two  input 
data  files:  DATAMAKE  DATA  and  HESSBLI  OUTPUT. 

DATAMAKE  DATA  is  the  number  of  points  on  the  body. 

Card  1  NN  (13) 

NN  Number  of  body  points. 

HESSBLI  OUTPUT  has  been  described  in  "Input  of  DATAMA." 

Description  of  Output  for  BLDMKER.  This  is  described  in  "Input 
to  Program  ICPL." 

Program  ICBLINT 

Description  of  the  Program  ICBLINT.  Program  ICBLINT  has  been 
developed  to 'provide  solutions  to  a  large  number  of  incompressible 
boundary-layer  flows.  The  program  can  calculate  either  body  only  or  body 
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and  near  wake  cases  and  can  treat  2-D  or  axi symmetric  geometries 
with  blunt  or  sharp  noses.  Solutions  may  be  obtained  using  only 
ICBLINT  or  ICBLINT  in  a  viscous-inviscid  iterative  procedure  with 
pressure  coefficients  of  the  displaced  body  as  the  medium  of  iteration. 
The  program  can  predict  laminar  and/or  turbulent  flows,  the  transition 
models  used  are  described  in  the  contractor  report  under  separate 
cover. 

The  flow  of  logic  in  the  program  is  controlled  by  MAIN  (see 
description  of  MAIN).  To  assist  the  user,  two  flow  diagrams  of  MAIN 
are  included.  Figure  8  gives  a  flow  chart  of  MAIN  by  the  logical 
functions  used,  while  Fig.  9  gives  a  flow  diagram  of  MAIN  by  the 
subroutines  called.  As  the  major  subroutines  have  particular  logical 
functions,  there  is  a  close  correlation  between  the  two  flow  charts. 

Additional  Comments  on  the  use  of  HESS/I CBLINT .  Although 
the  number  of  body  points  allowed  is  a  maximum  of  501,  for  the  body 
problems  of  concern  here,  fewer  than  101  points  are  preferred.  The 
reasons  for  this  are  twofold:  first,  only  a  maximum  101  profiles  may 
be  solved,  the  remainder  will  be  used  for  body  definition  only;  and 
second,  the  solutions  at  only  101  points  may  be  written  and  hence  used 
in  the  development  of  the  displaced  body  in  the  Potential  Flow  Code 
and  for  the  displaced  body  definition  in  ICBLINT.  For  separating 
flows,  the  number  of  body  points  should  be  checked  and. if  necessary, 
adjusted  after  CDCP  EXEC  has  been  executed. 

See  Appendix  A,  Table  3.5  for  a  listing  of  BLAYER  EXEC  and 
Table  3.6  for  a  listing  of  ICBLINT  FORTRAN. 
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Figure  8.  FLOW  DIAGRAM  OF  MAIN  ROUTINE  BY  FUNCTION 


START 


Figure  9.  FLOW  DIAGRAM  OF  MAIN  ROUTINE  BY 
SUBROUTINES  CALLED 


Description  of  Variables  in  Common.  In  this  section,  each 
variable  which  is  included  in  a  labeled  carmon  block  is  given  a  brief 


description.  The  name  of  the  labeled  common  block  in  which  the  variable 
occurs  is  given  to  the  left  of  the  variable  name. 


Common  Block 

Variable  Name 

Description  of  Variable 

NMLCRD 

ADTEST 

9 

Criterion  for  Increasing 
ETAINF,  (F( IE)  -  F( IE-4) 
is  compared  with  ADTEST) 

GEOME 

ANGLE 

9 

Local  Body  Angle,  in 
Degrees  (Converted  to 
Radians) 

TRANS 

ATR 

9 

Transition  Constant 
(=0.412) 

ARRAY  1 

AOB 

= 

Ao 

ARRAY 1 

AOBP 

ss 

3Aq 

3n 

ARRAY 1 

A1 

9 

Partial  Differential 
Equation  Coefficient 

COMWLL 

A1B 

9 

Derivative  of  U/Un  at 
Wall  e 

ARRAY 1 

A2 

9 

Partial  Differential 
Equation  Coefficient 

ARRAY 1 

A3 

9 

Partial  Differential 
Equation  Coefficient 

ARRAY 1 

A4 

9 

Partial  Differential 
Equation  Coefficient 

EDGPRP 

BETA 

S 

6  =  2c/Ue(dUe/dO 

CNVERG 

CCRNI 

« 

1  -  CRN I 

CFPR 

CF 

9 

Skin  Friction  in  Trans¬ 
formed  Coordinates 

(“  Cf  /EVD^ 
00 
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COEFF 

CFBAR 

= 

Cf  =  *  fo  Cf.  dt 

COEFF 

CFBREX 

= 

Cf  /REX 

COEFF 

CFE 

= 

C,  =  2t  /p  U2 
f  w'^e  e 

e 

COEFF 

CFINF 

= 

Cf  -  2t  Jp  U2 

T  W  00  00 

00 

COEFF 

CFRES 

= 

CFREY 

COEFF 

CFREY 

= 

cf  /Rex 

00 

COEFF3 

CF1 

9 

Term  used  in  Computing 
CFBREX 

ARRAY 1 

CHI 

= 

2 

x  ■  J-  fy  »  Vorticity 
Reynolds  Number 

TRANS 

CHI CRT 

9 

Value  of  CHIMAX  at  which 
Transition  Regime  Cal¬ 
culations  are  Begun 

TRANS 

CHIMAX 

* 

Maximum  Value  in  CHI 
Array 

CNVERG 

CONVRG 

* 

Convergence  Criterion 

CNVERG 

CRN  I 

9 

Selector  for  Finite 
Difference  Scheme;  1.0 
for  Fully  Implicit,  0.5 
for  Crank-Nicolson,  and 
0.0  for  Fully  Explicit 

COEFF 

DEL 

S 

6,  Boundary  Layer  Thick¬ 
ness  (Value  of  y  where 
U/Ue  =  0.995) 

COEFF2 

DELORF 

S 

6/reference  length 

COEFF2 

DELOX 

s 

6/x 

COEFF 

DELST 

s 

6*,  Boundary  Layer 
Displacement  Thickness 

COEFF 

DELTA 

s 

«/eVD 
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CNVERG 

DIF 

t 

Maximum  Difference  in  G 
or  F  (or  their  Derivative 
at  the  Wall)  between 
Iterations  (Compared  with 
CONVRG  in  Convergence 

Test) 

CNVERG 

DIFF 

> 

Local  Difference  in  G  or 

F  between  Iterations 

CNVERG 

DIF2 

» 

Differnece  in  the  Deriva¬ 
tive  of  F  at  the  Wall 
Between  Iterations 

ARRAY! 

DN 

> 

Array  of  An  Values 

ARRAY 1 

DN2 

» 

Array  of  An  Values  (Used 
only  when  ETAINF  is 
Increased) 

GEOME 

DS 

S 

AC  Corresponding  to  Ax 

COEFF2 

DSAXOR 

= 

6*  ./reference  length 

dX  1 

COEFF 

DSTARK 

= 

6£,  Incompressible  Bound¬ 
ary  Layer  Displace¬ 
ment  Thickness 

COEFF2 

DSTODL 

= 

6*/6 

COEFF2 

DSTORF 

= 

6* /reference  length 

COEFF2 

DSTOTH 

= 

<s*/e 

COEFF 2 

DSTOX 

= 

6*/x 

COEFF 

DSTRAX 

s 

6*  ,  =  r  {/l+26*/r  -  1} 

dX  i 

(Zero  for  2-D  Flow) 

EDGPRP 

DUEDS 

= 

dUE/ds 

GEOME 

DX 

= 

ax.  Step  Size  in  x 
Direction 

GEOME 

DXMAX 

> 

Maximum  Value  of  ax  Per¬ 
mitted  in  the  Calculations 

GEOME 

DXOLD 

» 

Previous  Value  of  ax 

GEOME 

DX! 

» 

Storage  Value  of  ax  at 
the  First  Station 

GEOME 


DYDX 


GEOME 

DYDX 

• 

Body  Slope  at  JSTA 

ARRAY! 

EPSO 

= 

£q,  the  Outer  Eddy  Viscosity 

ARRAY1 

EPSPL 

= 

et,  or  E+,  the  Inner  Eddy 

Viscosity  or  the  Combined 
Inner  and  Outer  Eddy 
Viscosities 

VSCSTY 

EPSVD 

= 

eVD  =  {uref/(p-U“L)}  2’  where 

L  is  the  Unit  Length  in  REINF 

NMLCRD 

ETAINF 

• 

Maximum  Value  of  n 

COMWLL 

El 

• 

Matrix  Inversion  Coef¬ 
ficient  Evaluated  at  the  Wall 

ARRAY 1 

FC 

= 

CRN I ( F2-F1 )  +  FI 

ARRAY! 

FCN 

= 

CRN  I ( F2N -FIN )  +  FIN 

ARRAY! 

FCP 

= 

(F2-F1)/DS 

COMWLL 

FF 

» 

Wall  Boundary  Condition 
on  F 

ARRAY 1 

FI 

* 

Value  of  F  (U/Ug )  from 

Previous  Iteration  or  Station 

ARRAY 1 

FIN 

= 

aFl 

3n 

ARRAY 1 

FINN 

__ 

32F1 

2 

3n^ 

ARRAY 1 

F2 

9 

Value  of  F  from  Solution 
of  Momentum  Equation 

ARRAY 1 

F2N 

= 

aF2 

3n 

ARRAY 1 

F2NN 

32F2 

3n2 

COMWLL 

F2N1 

• 

Value  of  F2N  at  Wall 
from  Previous  Iteration 

27 


TRANS 

GAMMA 

9 

Transition  Intermittency 
Factor,  r,  r  =  1  -  exp 
—2 

( -As  ) ,  where  A  =  ATR  and 
?  =  XI BAR 

COEFF 

HAFCF 

= 

Cf  /2 

00 

INTGR 

IE 

9 

Number  of  Grid  Points  in 
n  Direction 

NTEGER 

II 

9 

Subscript  for  Array  IPR 

INTGR 

UMAX 

9 

Maximum  Number  of  Values 
in  Geometry  Arrays  (i.e. 

ZA,  XSTA,  etc.) 

INTGR 

IM 

= 

IE  -  1 

NTEGER 

IPFL 

9 

Number  of  Values  in  Array 
IPRFL 

ARRAY2 

IPR 

9 

Array  of  Subscripts  for 

Array  XSTA  Giving  Values 
of  x  at  which  a  Solution 
is  to  be  Obtained 

ARRAY2 

IPRFL 

9 

Array  of  Subscripts  for 

Array  XSTA  Giving  Values 
of  x  at  which  Complete 
Solutions  ( with  Profiles) 
are  to  be  Printed 

NTEGER 

I  PRINT 

9 

Number  of  Values  in  Array 

IPR 

INTGR 

IPRNT 

9 

Subscript  for  Array  IPRFL 

INTGR 

I  STOP 

9 

Number  of  Iterations  Since 
Last  Converged  Solution 

NTEGER 

JJ 

9 

Integer  Used  in  Computing  DX 

APLDAT 

JSTA 

9 

Beginning  of  N-S  Box 

NTEGER 

K 

9 

Station  Counter 

NTEGER 

KADETA 

9 

Indicator,  if  Zero,  ETAINF 
is  Held  Constant,  if 
Non-zero,  ETAINF  can  be 
Increased 
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NTEGER 

KEND 

9 

Maximum  Value  Allowed 
for  K 

NTEGER 

KEP 

9 

Subscript  in  EPSPL  Array 
Where  the  Inner  and  Outer 
Eddy  Viscosities  are 
Matched 

INTGR 

KL 

9 

Index  for  Profile  Output 

DO  Loop,  1  Gives  Output 
for  Each  Grid  Point,  2 
for  Every  Second  Grid 

Point 

NTEGER 

KSTOP 

9 

Indicator,  Normally  Zero, 
Set  to  One  when  Solution 
is  Completed  (X  =  XSTA 
(UMAX)) 

TRANS 

KTRANS 

9 

Indicator  for  Transition 
Regime,  Zero  for  No 
Transition  or  Instantane¬ 
ous  Transition,  One  for 
Transition  Regime  (Set  at 
Two  at  "End  of  Transition" 

NTEGER 

KTRNSN 

9 

Subscript  in  XSTA  Array 
at  Which  Instantaneous 
Transition  Occurs 

INTGR 

KVSLAW 

9 

Indicator  for  Inner  Eddy 
Viscosity  Law,  Zero  for 
Reichardt  Law,  Non-zero 
for  Van  Driest  Law 

NTEGER 

LAMTRB 

9 

Indicator,  1  for  Laminar 
Flow,  2  for  Turbulent 

Flow 

CNVERG 

NC 

9 

Specifies  Convergence  Test 
Zero  for  Test  on  Deriva¬ 
tives  at  Wall ,  One  for 

Test  on  Functions  at  all 
Points  in  Profile 

INTGR 

NIT 

9 

Number  of  Iterations  at 
Current  Value  of  x 

NTEGER 

NITTOT 

9 

Cumulative  Total  Number 
of  Iterations  in  x 
Direction 
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INTGR 

NIT1 

9 

If,  after  a  converged 
solution  is  obtained, 

NIT  is  less  than  or  equal 
to  NIT1,  ax  is  doubled 

INTGR 

NIT2 

9 

If,  after  a  converged 
solution  is  obtained,  NIT 
is  greater  than  or  equal 
to  NIT2,  ax  is  halved 

INTGR 

NIT3 

9 

If  NIT  =  NIT3  before  a 
converged  solution  is 
obtained,  ax  is  halved, 

X  and  XI  are  recomputed 
and  the  iteration  loop  is 
restarted  with  the  smaller 
values  of  X  and  XI 

APLDAT 

NOPRFS 

9 

Number  of  Profiles  in 

N-S  Box 

INTGR 

NOSE 

9 

Indicator  of  Geometry,  1 
for  a  Blunt  Nose,  2  for 
a  Sharp  Nose 

STRT 

NRSTRT 

9 

Increment  in  K  for  which 
the  Restart  Tape  is 

Written 

EDGPRP 

PE 

= 

Vref-  "hera 

Pref  "  pref  ^ref 

EDGPRP 

PESO 

= 

Pe/Pref-  at  5  *  0 

FRSTRM 

PFS 

9 

Freestream  Static 
Pressure,  PS  A 

GEOME 

PNC 

/?7 

— 12 — j  ,  Term  in  Trans- 

u  r^ 
e  w 

formation  of  Normal 
Coordinate 

EDGPRP 

PP 

= 

dPE/ds 

REF 

PREF 

= 

2 

^ref  ~  pref  ^ref 

STAG 

P  ST  AG 

9 

Freestream  Stagnation 
Pressure,  PSFA 

ARRAY2 

PZ 

* 

Local  Pressure  Array,  P/Prpf 

STAG 

PIO 

= 

VPref 

GEOME 

REFLEN 

9 

Reference  Length  (may  be 
any  Value  Appropriate  to 
User's  Needs 

FRSTRM 

REINF 

9 

Freestream  Reynolds  Num¬ 
ber  per  Unit  Length 

COEFF 

RETHET 

= 

Re0,  Edge  Unit  Reynolds 

Number  Times  e 

COEFF 

REX 

9 

Local  Reynolds  Number 

Based  on  Edge  Conditions 

FRSTRM 

RHOFS 

= 

p  ,  slugs/ft3 

REF 

RHOREF 

= 

pref  =  Pc’  slugs/ft3 

C0EFF2 

ROREFL 

= 

r/reference  length 

ARRAY2 

RZ 

= 

r.  Array  of  Radius  Values 

GEOME 

RO 

= 

r.  Local  Radius 

GEOME 

SCF 

9 

Scale  Factor  used  to  Con¬ 
vert  Unit  Length  for  REINF 
(REINF  is  Computed  per 

Foot  if  SCF=1;  per  inch 
if  SCF=12 . ) 

C0EFF3 

SUM 

9 

Term  used  in  Computing 

CFBREX 

COEFF 

THET 

= 

e.  Boundary  Layer  Momen¬ 
tum  Thickness 

COEFF2 

THODEL 

= 

0/6 

COEFF2 

THOREF 

= 

©/Reference  Length 

EDGPRP 

UE 

= 

Wef 

EDGPRP 

UER02 

= 

<ue/uref>r2J 

EDGPRP 

UESO 

S 

Ue/Uref  at  s  =  0 
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ARRAY2 

UEZ 

= 

Ue/U-’  Arrdy 

FRSTRM 

UFS 

= 

U  ,  ft/ sec 

00 

REF 

UREF 

= 

uref  =  U«’  ft/sec 

ARRAY 1 

VC 

= 

'•-£  <-2s  tf -F><in. 
Transformed  Normal 

Velocity  Component 

COMWLL 

VW 

9 

Value  of  VC  at  Wall 

GEOME 

X 

9 

Local  Surface  Distance,  x 

TRANS 

XBAR 

t 

9 

Measure  of  Relative  Length 
of  Transition 

^*e"  XEnd  of  Transition" 

XBAR* *Beg inning  of  Transi¬ 
tion) 

GEOME 

XI 

= 

£,  Transformed  Surface 
Distance 

TRANS 

XI  BAR 

= 

X,  Transition  Inter- 

mittency  Factor  Coordinate 
X  =  (X-XZERO)/XLAMDA, 
where  XLAMDA  =  (XBAR-1)* 
XZERO/4  and  XZERO  is  the 
Value  of  x  at  the  Start  of 
Transition 


GEOME 

XI  OLD 

9 

Previous  Value  of  XI 

GEOME 

X12 

= 

2*XI 

GEOME 

XJAY 

= 

XJFAC 

GEOME 

XJFAC 

9 

Indicator,  Zero  for  2-D 

Flow,  One  for  Axially 
Symmetric  Flow 

NMLCRD 

XKETA 

9 

Parameter  Which  Controls  the 

Grid  Spacing  in  the  n  Direc¬ 
tion  (1.0  Gives  an  Equally 
Spaced  Grid,  a  Value  Greater 
than  1.0  Gives  a  Grid  with 
a  Smaller  Spacing  at  the  Wall 
than  at  the  Outer  Edge) 


VSCSTY 

XK1 

9 

Constant  for  Van  Driest 
Inner  Eddy  Viscosity  Law 

VSCSTY 

XK2 

9 

Constant  for  Outer  Eddy 
Viscosity  Law 

ARRAY1 

XN 

9 

Array  of  n  Values 

ARRAY1 

XN2 

9 

Array  of  n  Values  (Used 
only  when  ETAINF  is  In¬ 
creased) 

GEOME 

XOLD 

9 

Previous  Value  of  X 

C0EFF2 

XOREFL 

= 

x/Reference  Length 

ARRAY2 

XSTA 

= 

x.  Array  of  Surface  Dis¬ 
tance  Values  used  to 
Specify  Geometry 

APRAY1 

Y 

S 

y/£vD.  Stretched  Normal 
Coordinates 

ARRAY 1 

YOVDEL 

= 

y/6 

ARRAY1 

YOVTHT 

= 

y/e 

ARRAY 1 

YY 

= 

y.  Normal  Coordinate 

GEOME 

Z 

= 

z.  Local  Axial  Coordinate 

ARRAY2 

ZA 

= 

z.  Array  of  Axial  Co¬ 
ordinate  Values  used  to 
Specify  Geometry 

GEOME 

ZOL 

= 

z/Body  Length 

COEFF2 

ZOREFL 

= 

z/Reference  Length 

COEFF3 

SUM1 

Incremental  Viscous  Drag 
due  to  Finite  Surface 
Distance 

COEFF3 

SUMT 

Total  Drag  due  to  Skin 
Friction  and  Pressure 

STAG 

RHO 

Density  from  Stagnation/ 
Static  Pressure 

TVCURV 


TVC 


Flag  for  TVC  Calculations 
Yes  for  TVC,  No  for  2-D 
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TVCURV 

BODANG 

Local  Body  Angle 

TVCURV 

RTVC 

Local  Body  Radius 

TVCURV 

RORW 

TVC  Body  Radius 

TVCURV 

RORWSQ 

RORW 

TVCURV 

RW 

Local  Body  Radius,  ft. 

SPTVC 

RNOSE 

Nose  Radius,  ft. 

WTE 

IW 

Counter  for  IWT 

WTE 

IWT 

Increment  for  Printing 

Full  Profile  Data 

TE 

KTE 

Station  Number  of  '.'rail¬ 
ing  Edge 

TE 

KTSTTE 

Flag  for  L.E.  Station 
Calculations 

WAK 

XKFAC 

Flag  for  Making  Wake 
Calculations;  0  for  No  Wake 
1  Include  Wake 

HESS 

IWRTH 

Displacement  Body-Original 
Body  Weighting  Factor 

HESS 

ICARD1 

Hess  Code  Execution 
Parameter,  Printed  on 
Bodycoor  Data 

HESS 

ICARD2 

HESS 

ICARD3 

HESS 

ICARD4 

HESS 

NN 

Number  of  Body  Stations 

HESS 

IBDN 

Body  Sequence 

HESS 

HEDR 

Array  of  Alphanumeric 
Characters  Describing 
the  Executing  Case 

HESS 

CASE 

Alphanumeric  for  Record¬ 
ing  Case  Number 
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HESS 

KHESS 

Flag:  1  if  Inviscid  C 

Data  are  Read  in  from  p 

HESS,  0  if  no  Inviscid 

Cp  Data  are  Read  in 

DRAG 

CDV1 

Old  Value  of  CDV2 

DRAG 

CDV2 

Pressure  Drag  for  In¬ 
cremental  Length  dx 

DRAG 

XCP 

Cp  at  x  from  Current 

Iteration 

DRAG 

XCP1 

Cp  at  x  from  Previous 

Iteration 

DRAG 

CP 

Cp  Array  from  Current 

Iteration 

DRAG 

CPI 

Cp  Array  from  Previous 
Iteration 

DRAG 

ROMAX 

Maximum  Body  Radius 

INTACT 

BETAS 

Minimum  BETA 

STRTUP 

I  READ 

Flag  for  Restart  Procedure; 

1  use  Restart,  0  do  not  use 
Restart 

STRTUP 

ITER 

Station  Number  at  Which 

Restart  Tape  is  Written 

STRTUP 

ITER 

Iteration  Number 

BDCOOR 

ZDELST 

Axial  Displaced  Body  Co¬ 
ordinate  at  Original  Axial 
Coordinate  Station 

BDCOOR 

RODELS 

Displaced  Body  Radius  at 
Station  k 

BDCOOR 

ZOLD 

Value  of  Z  from  Previous 
Iteration 

BDCOOR 

ROLD 

Value  of  R  from  Previous 
Iteration 

BDCOOR 

MM 

IRST 

APLDAT 

XLOCAT 

Axial  Location  of  the  Upstream 
N-S  Region  Boundary 
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APLDAT 

APLDAT 

APLDAT 


DYDX 

JSTA 

NOPRFS 


Body  Slope  over  the  Tail 
Section,  Taken  as  Constant 
for  the  Straight  Taper  Re¬ 
quirement  of  the  APLNS 
Code 

Station  Number  of  the  Be¬ 
ginning  of  the  N-S  Box 

Number  of  APLNS  Code 
Stations  on  the  Body 


Description  of  Subroutines. 

Subroutine  MAIN.  The  principle  function  of  MAIN  is  to 
control  the  flow  of  logic  in  the  program-,  however,  some  calculations  are 
also  made  in  MAIN.  In  general,  the  calculations  in  MAIN  are  such  that 
they  cannot  conveniently  be  performed  in  the  subroutines  or  they  are 
not  sufficiently  similar  to  the  calculations  made  in  particular  sub¬ 
routines. 

Subroutine  READIN  is  called  first.  If  the  input  data  are  from 
the  restart  tape,  the  calculations  are  made  which  normally  occur  at  the 
end  of  MAIN  after  subroutine  WRITE2  is  called,  i.e.,  the  calculations  that 
enable  the  program  to  set  up  for  the  next  station  are  made.  The  next 
section  of  the  program  is  skipped,  and  the  marching  integration  is  re¬ 
sumed.  Otherwise,  the  initial  data  and  profiles  are  calculated  and 
written  by  subroutine  WRITE1  and  the  integration  is  begun  with  the  state- 
ment;  DO  330  K=KK,  KEND.  An  iteration  loop  is  then  entered  for  the  solu¬ 
tion  of  the  velocity  profile  at  that  station.  The  boundary  conditions 
at  the  wall  are  calculated  and  the  momentum  equation  is  solved.  The 
transformed  normal  velocity  component  V  is  calculated  together  with  the 
skin  friction  in  the  transformed  coordinates  (CF).  The  iteration  counters 
are  incremented  and  NIT  and  ISTOP  are  checked.  If  ISTOP  is  greater  than 
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100,  the  program  is  stopped.  If  NIT  equals  (or  exceeds)  NIT3  the  step- 
size  is  halved,  X  and  DX  are  recalculated,  and  control  is  returned  to 
the  beginning  of  the  iteration  loop.  Otherwise,  the  next  test  is  the 
convergence  test.  If  the  convergence  criterion  is  not  satisfied,  con¬ 
trol  is  returned  to  the  beginning  of  the  iteration  loop. 

When  convergence  has  been  obtained,  the  velocity  profile  is 
checked  to  see  if  n  needs  to  be  increased.  If  so,  subroutine  ADDETA 

oo 

is  called  and  control  is  returned  to  the  beginning  of  the  iteration  loop. 
If  not,  subroutines  COEF  and  WRITE2  are  called.  New  values  of  x  and  s 
are  obtained  and  the  present  values  of  velocity  are  assigned  to  the 
arrays  which  contain  the  values  from  the  previous  iteration,  and  the  it¬ 
eration  counters  NIT  and  ISTOP  are  reset.  If  a  solution  has  been  ob¬ 
tained  for  the  final  value  of  x,  (or  for  K  =  KEND)  the  program  stops, 
otherwise,  the  program  returns  to  the  beginning  of  the  marching  inte¬ 
gration  loop  to  obtain  a  solution  at  the  new  value  of  x. 

Description  of  Variables. 

ABST  =  F2 ( IE )  -  F2( IE-4) .  If  ABST  is  greater  than  or  equal  to 

ADTEST,  n  is  increased 

co 

IPRIIM  =  IPR  (II-l) 

KREADN  ;  Indicator  for  Subroutine  READIN,  1  for  Input  Data,  2  for 
Initial  Data  Written  on  Restart  Tape,  3  for  Station  Data 
Written  on  Restart  Tape 
N  ,  Array  subscript 

The  following  variables  used  in  MAIN  are  described  in  the  section 
Description  of  Variables  in  Common:  ADTEST,  A1B,  CCRNI,  CF,  CF1,  CONVERG 
CRN I ,  DIF,  DN,  DUEDS,  DX,  EPSVD,  FC,  FCN,  FCP,  FI,  FIN,  F2,  F2N,  F2NN, 
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F2N1  IE,  II,  UMAX,  IPR,  IPRNT,  ISTOP,  K,  KADETA,  KEND,  KSTOP,  KSTRT, 
LAMTRB,  NIT,  NITTOT,  NIT3,  NOSE,  NT3,  PNC,  REINF,  RHOFS,  RO,  SCF,  UE, 

UFS,  UREF,  VC  VW,  X,  XI,  XIOLD,  XI2,  XJFAC,  XKFAC,  XN,  XOLD,  XSTA,  Y, 

YY. 

Subroutine  ADDETA.  Subroutine  ADDETA  is  used  to  increase/ 
decrease  ETAINF  when  the  value  of  F  at  the  outermost  grid  point  and  the 
value  of  F  four  points  inward  differ  by  more/less  than  the  quantity  ADTEST. 
A  new  array  of  values  is  generated  with  the  maximum  value  increased  by 
10  percent.  Three  point  interpolation  is  used  to  obtain  the  values  of 
velocity  at  the  new  grid  points,  whereas  at  the  grid  points  lying  between 
the  new  and  old  values  of  f  is  set  equal  to  unity.  First  and  second 
derivatives  of  F  are  then  calculated.  The  storage  arrays  for  F  and  their 
derivatives  are  updated  and  V  is  recalculated  using  the  new  values.  The 
arrays  for  y  are  not  updated  in  this  subroutine  but  are  recalculated  at 
the  normal  points  in  the  solution  procedure. 

The  subroutine  provides  the  output  of  the  intermediate  profiles 
following  the  message  that  n,,,  has  been  increased.  DX,  XKETA,  ETAINF, 
and  ADTEST  are  also  included  in  the  output. 

USAGE: 

CALL  ADDETA 


Description  of  Variables. 

DETA1 

An  at  wall 

ETAIN2 

New  value  n 

00 

JC 

Array  subscript 

N 

Array  subscript 

38 


The  following  variables  used  in  subroutine  ADDETA  are  described 
in  the  section  Description  of  Variables  in  Conmon:  AOBP,  CHI,  DN,  DN2, 

DX,  EPSPL,  ETAINF,  FC,  FCN,  FCP,  FI,  FIN,  FINN,  F2,  F2N,  F2NN,  IE,  IM, 
ISTOP  Kl,  VC,  XKETA,  XN,  XN2,  YOVDEL,  YOVTHT,  YY. 

Subroutine  BLUNT1 .  This  subroutine  calculates  the  edge 
and  reference  conditions  at  the  stagnation  point  for  a  blunt  body.  The 
derivative  of  Ue  with  respect  to  s  is  obtained  using  the  five-point  dif¬ 
ferentiation  formula  with  two  points  reflected  about  the  stagnation  point. 
USAGE: 

CALL  BLUNT 1 

Description  of  Variables. 

DUM  Dummy  variable 

DUM1  Dummy  variable 

DUM2  Dummy  variable 

DUM3  Dummy  variable 

J  Array  subscript 

NERR  Dummy  variable 

PE  Array  of  PL  values 

UE  Array  of  U  values 

The  following  variables  used  in  Subroutine  BLUNT1  are  described 
in  the  section  Description  of  Variables  in  Common:  BETA,  DUEDS,  PE, 

PESO  PFS,  PNC,  PP,  PREF,  PSTAG,  PZ,  POPRIM,  P10,  REINF,  RO,  SCF,  UE, 

UER02  UFS,  UREF,  X,  XJFAC,  XM,  XSTA,  Z,  ZA. 

Subroutine  BLUNT2.  For  a  blunt  body,  the  edge  properties 
required  by  subroutine  DELTAS  to  calculate  AC  are  provided  by  subroutine 
BLUNT2.  In  the  immediate  vicinity  of  the  stagnation  point,  z,  r,  Pe  and 


dUe/ds  are  calculated  using  five  point  interpolation  with  two  points 
reflected  about  the  stagnation  point. 

USAGE: 

CALL  BLUNT2 

Description  of  Variables. 


DUM 

Dummy  variable 

DUM1 

Dummy  variable 

DUM3 

Dummy  variable 

J 

Array  subscript 

The  following  variables  used  in  subroutine  BLUNT2  are  described 
in  the  section  Description  of  Variables  in  Common:  DUEDS,  IIAX,  PE,  PP 
PZ,  P10,  RZ,  RO,  UE,  UEZ,  X,  XSTA,  Z,  ZA. 

Subroutine  CHANGE.  Subroutine  CHANGE  adjusts  the  x  step- 
size  when  necessary  to  obtain  a  solution  at  the  specified  values  of  x. 
Subroutine  CHANGE  is  also  used  when  calculations  are  made  with  the  in¬ 
stantaneous  transition  model.  In  this  case  LAMTRB  is  reset  to  two,  the 
stepsize  is  reduced  and  a  message  is  printed. 

Subroutine  CHANGE  also  increases  or  decreases  the  x-stepsize 
depending  on  the  number  of  iterations  required  to  converge  the  previous 
station. 

USAGE: 

CALL  CHANGE 

Description  of  Variables. 

IPRII  =  I  PR ( 1 1 ) 

IPRIIM  =  IPR(II-l) 

The  following  variables  used  in  subroutine  CHANGE  are  described 
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in  the  section  Description  of  Variables  in  Common:  DX,  DXMAX,  DXOLD, 

DX1 ,  EPSVD, ,  II,  UMAX,  IPR,  JJ ,  KSTOP,  DTRNSN,  LAMTRB,  NIT,  NIT1, 

NIT2,  X,  XOLD,  XSTA. 

Subroutine  COEF.  Subroutine  COEF  calculates  the  skin- 
friction  coefficients  for  the  boundary-layer  flow,  the  pressure  drag, 
and  total  drag.  It  calculates  the  Reynolds  number,  displacement  thick¬ 
ness,  momentum  thickness  and  other  dimensionless  parameters.  At  X  =  0 
the  calculations  in  this  subroutine  are  omitted  except  for  Z/L. 

USAGE: 

CALL  COEF 

Description  of  Variables. 

CF2  Term  used  in  computing  CFBREX 

DSUM  Term  used  in  computing  CFBREX 

FAC1  Dummy  variables 

FAC2  Dummy  variables 

I  Array  subscript 

N  Array  subscript 

RESE  Re„4 

SQREX  Square  root  of  REX 

STRES  CHREY 

The  following  variables  used  in  subroutine  COEF  are  described  in 
the  section  Description  of  Variables  in  Common:  A1B,  CF,  CFBAR,  CFBREX, 
CFINF,  CFRES,  CFREY,  CF1,  DEL,  DELORF,  DELOX,  DELST,  DELTA,  DN,  DSAXOR, 
DSTODL,  DSTORF,  DSTOTH,  DSTOX,  STRAX,  DX,  EPSVD,  FC,  HAFCF,  IE,  UMAX  K, 
KFS,  LAMTRB,  REFLEN,  RE  INF,  RETHET,  REX,  PNC,  RHOFS,  ROREFL,  RO  SUM, 
THET,  THODEL,  THOREF,  UE,  UFS,  X,  XJFAC,  XOREFL,  Y,  YOVDEL,  YOVTHT,  YY, 


Z,  ZA,  ZOL,  ZOREFL. 

Subroutine  DELTAS.  Subroutine  DELTAS  calculates  the  trans¬ 
formed  surface  distance  £  using  integration  by  Simpson’s  rule.  A  value 
of  £  is  obtained  for  each  increment  A£  from  the  equation 

d£  =  U  r2Vkdx 
ewe 

The  transformation  factor  PNC, 


rnv,  ~  - t- 

U  rJ 
e 

and  the  pressure  gradient  term,  $  , 

-  21 
Ue  d£ 

are  also  calculated. 

USAGE: 

CALL  DELTAS 

Description  of  Variables. 


HA 

X  -  DX/2 

RP 

RO  at  X 

ROHA 

UE  at  HA 

XI 

Storage  value  of  X 

The  following  variables  used  in  subroutine  DELTAS  are  described 
in  the  section  Description  of  Variables  in  Common:  BETA,  DS,  DUEDS,  DX, 
PNC,  RO,  UE,  UER02,  X,  XI,  XIOLD,  X12,  XJFAC,  XKFAC. 

Subroutine  DERIV3.  Subroutine  DERIV3  uses  subroutine  FD3 


to  generate  the  first  derivative  of  a  function  which  is  given  by  the 


array  F.  The  array  of  abscissa  values  is  given  by  X  and  the  derivative  of 
F  with  respect  to  X  is  returned  in  the  array  FP. 

USAGE: 

CALL  DERIV3  (F,X,IMAX,IMIN,FP)  where  F  and  X  are  the  arrays  of 

ordinate  and  abscissa  values,  and  where  IMAX  and  IMIN  are  the  upper  and 

lower  subscripts  for  the  array  FP. 

Subroutine  EFFMU.  Subroutine  EFFMU  calculates  the  terms 

in  the  coefficients  of  the  differential  equations  which  contain  the  eddy 

3Ag 

viscosity.  These  terms  are  Aq  and  where: 

Aq  =  (l+e+) 

for  the  momentum  equation.  For  the  laminar  case,  c+  =  0  and  the  terms 
reduce  to 

Aq  =  1 
and 

Aq  "  0 

The  first  portion  of  the  subroutine  calculates  the  terms  for  the 
laminar  case,  and  the  next  portion  of  the  subroutine  calculates  the 
vorticity  Reynolds  number,  x, 

2 

X=L  i y. 

V  3y 

For  calculations  with  a  transition  regime,  the  routine  initiates 
transition  and  calculates  the  transition  intermittency  factor  r  where 
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and  Xg  is  the  location  of  the  beginning  of  transition  and  x  is  a  measure 
of  the  relative  length  of  transition. 

In  the  third  portion  of  the  subroutine,  the  eddy  viscosity  is 

l 

computed  for  turbulent  flow  calculations  and  the  terms  Ag  and  Ag  are 
then  computed. 

USAGE: 

CALL  EFFMU  (LAMTRB,  II,  DSTARK) 

Description  of  Variables. 

ARG  Dummy  variable 

DELTA  6/eyg,  boundary-layer  thickness  in  stretched  coordinates 

„  ★ 

DSTARK  <$k,  incompressible  displacement  thickness 

FAC  Dummy  variable 

FAC1  Dummy  variable 

FAC2  Dummy  variable 

FAC3  Dummy  variable 

FAC4  Dummy  variable 

GAMMA  y,  outer  eddy  viscosity  intermittency  factor 

I  The  subscript  in  the  EPSPL  array  where  the  inner  and  outer 

eddy  viscosities  are  matched. 

LAMTRB  Indicator  for  laminar  or  turbulent  flow;  1  for  laminar, 

2  for  turbulent 

N  Array  subscript 

XLAMDA  X  =  (x-l)x0/4 

XZERO  Xg,  location  of  start  of  transition 

The  following  variables  used  in  subroutine  EFFMU  are  described 
in  the  section  Description  of  Variables  in  Common:  ATR,  AOB,  AOBP,  CF, 
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CHI,  CHI CRT,  CHIMAX,  DN,  DX,  DXOLD,  EPSO,  EPSPL,  EPSVD,  FC,  FCN,  GAMMA, 
PRT,  RO,  UE,  X,  XBAR,  XI BAR,  XI2,  XJFAC,  XKFAC,  XK1 ,  XK2 ,  XN,  Y. 

Subroutine  EDGPROP.  Subroutine  EDGPROP  obtains  the  edge 
properties  used  in  subroutine  DELTAS  to  calculate  A£.  The  edge  properties 
are  obtained  from  subroutine  BLUNT1. 

USAGE: 

CALL  EDGPROP 

Subroutine  FD3  and  FD5.  Subroutines  FD3  and  FD5  use  the 
first  derivative  of  the  Lagrangian  interpolating  polynomial  of  second 
and  fourth  order,  respectively,  to  provide  the  first  derivative  of  the 
function  F  at  point  X.  The  returned  value  is  denoted  by  FX.  The  inter¬ 
polating  polynomials  are  described  in  the  descriptions  of  subroutines 
INTER3  and  INTER5. 

USAGE: 

CALL  FD3  (X,X1 ,X2,X3,F1 ,F2,F3,FX)  or 

CALL  FD5  (X,X1  ,X2,X3,X4,X5,F1 ,F2,F3,F4,F5,FX)  where  XI,  X2,  etc. 
are  the  abscissa  values,  and  where  FI,  F2,  etc.  are  the  ordinate  values. 

Subroutine  GEOM.  This  subroutine  obtains  the  edge 
properties  at  the  initial  value  of  x  for  the  geometry  specified  by  the 
variable  NOSE  which  is  passed  through  common.  The  edge  properties  are 
obtained  from  subroutine  BLUNT1  for  a  blunt  nosed  body  (NOSE  =1), 
from  subroutine  CONE!  for  a  sharp  nose  (NOSE  =2). 

USAGE: 

CALL  GEOM 

Subroutine  I N IT .  Many  variables  are  given  fixed  initial 


values  in  the  program  before  the  marching  integration  procedure  is  begun. 
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The  non-subscripted  variables  are  given  the  initial  values  in  subroutine 
INIT. 

USAGE: 

CALL  INIT 

Description  of  Variables. 

JT  Array  subscript 

The  following  variables  used  in  subroutine  INIT  are  described  in 
the  section  Description  of  Variables  in  Common:  CCRNI,  CF,  CHIMAX,  CRNI, 
DS,  DX,  DXOLD,  DX1 ,  El,  FF,  GAMMA,  IE,  UMAX,  IM,  IPRNT,  ISTOP,  JJ , 

KEP,  KFS,  KSTOP,  KTPW,  KTW,  LAMTRB ,  NIT,  NITTOT,  NOSE,  NRSTRT ,  RHOFS, 

SUM,  UFS,  VW,  X,  XI,  XI BAR,  XIOLD,  XI2,  XJAY,  XJFAC,  XKFAC,  XM,  XOLD, 
XSTA,  Z. 

Subroutine  INTER3.  Subroutine  INTER3  uses  a  second-order 
Lagrangian  interpolating  polynomial,  interpolating  on  the  points  xl,  x2, 
and  x3,  with  the  corresponding  function  values  FI,  F2,  F3  to  provide  a 
value  F(x). 

The  general  form  of  the  polynomial  is 
f 

where 

3 

m*)  -  n 

m=  1 
m^k 

USAGE: 

CALL  INTER3  (X,  XI ,X2,X3,F1 ,F2,F3,F)  where  XI,  X2,  and  X3  are 
the  abscissa  values  and  FI,  F2,  and  F3  are  the  function  values. 


A 
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Subroutine  INTERS.  Subroutine  INTER5  uses  a  fourth-order 
Lagrangian  interpolating  polynomial  interpolating  on  the  points  xl,  x2, 
x3,  x4,  and  x5  with  the  corresponding  function  values  FI,  F2,  F3,  F4, 
and  F5  to  provide  a  value  F(x). 

The  general  form  of  the  polynomial  is 

p(x)  *  X^]  /xLk(x) 

k=l  J  x  K 


where 


Lk(x)  = 


USAGE: 


CALL  INTER5  (X,X1 ,X2,X3,X4,X5,F1 ,F2,F3,F4,F5,F)  where  XI,  X2,  etc. 
are  the  abscissa  values  and  FI,  F2,  etc.  are  the  function  values. 

Subroutine  INTERP.  Subroutine  INTERP  uses  function  TLU 
to  linearly  interpolate  in  the  array  F2  for  the  value  FF  corresponding 
to  the  value  XX  in  the  array  XN.  If  XX  is  not  within  the  range  of  XN, 
a  message  is  printed  and  FF  is  set  equal  to  F2(IE)  where  IE  is  the  size 
of  the  arrays.  XN  must  be  a  strictly  increasing  array. 

USAGE: 

CALL  INTERP  ( XX , XN , F2 , 1 E , FF ) 

Subroutine  INTERP5.  Subroutine  INTRP5  uses  five-point 
Interpolation  to  obtain  a  value  FO  in  the  array  FI  corresponding  to  the 
value  XO  in  the  array  XI.  The  array  XI  can  be  either  monotone  increasing 
or  monotone  decreasing. 


i 
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USAGE: 

CALL  INTRP5  (XO, XI , FI , IMAX.IMIN.il ,F0)  where  XI  and  FI  are  the 
arrays  of  coordinate  values,  IMAX  and  IMIN  are  the  maximum  and  minimum 
values  allowed  for  the  subscripts  of  the  arrays  0  and  II  is  an  indicator 
(II <0 ,  XI  is  monotone  decreasing;  1 1 >0 ,  XI  is  monotone  increasing). 

Subroutine  MOMENT.  The  solution  of  the  momentum  equation 
is  obtained  in  subroutine  MOMENT.  The  coefficients  of  the  partial  dif¬ 
ferential  equation  are  calculated  and  new  values  of  F  =  U/Ug  and  the 
first  and  second  derivatives  of  F  with  respect  to  n  are  obtained  from 
subroutine  SOLVE.  The  difference  between  the  new  and  the  former  values 
is  obtained  for  either  the  wall  or  the  all  points  convergence  test.  The 
derivative  of  F  with  respect  to  x  is  also  calculated. 

USAGE: 

CALL  MOMENT 

The  following  variables  in  subroutine  MOMENT  are  described  in  the 
section  Description  of  Variables  in  Common:  AO,  AOB,  AOBP,  A1 ,  A 2,  A3, 
A4,  BETA,  C,  CCRNI ,  CP;  CRNI,  DIFF,  DS,  DSTARK,  El,  FC,  FCN,  FCP,  FF, 

FI,  FIN,  FINN,  F2 ,  F2NN,  F2N1,  IE,  IM,  KEP,  LAMTRB,  NC,  TH,  VC,  X,  XI. 

Subroutine  PRFILE.  Subroutine  PRFILE  calculates  the  array 
of  n  values,  XN,  which  corresponds  to  the  values  of  IE,  ETAINF,  and 

XKETA.  The  grid  spacing  is  given  by  Ani  =  KAn^  where  K  corresponds  to 

XKETA  and  thus  at  the  ith  grid  point 

ni  =  Anl  £ IT  *  1  =  °*  2*  3’  •••N  (N  =  IE  -  1) 

where  An-j  is  given  by 

a  -  K-l 

ni  "  n~  7J 
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The  initial  velocity  profile,  F,  is  calculated  from  the  formula 
F  =  1  -  exp(-n).  The  derivatives  of  F  are  also  calculated  as  are  initial 
profiles  for  e+,  y/5,  y/e.  Further,  the  boundary  conditions  at  the  outer 
edge  of  the  boundary  layer  are  set  in  this  subroutine. 

USAGE: 

CALL  PRFILE 

Description  of  Variables. 

DETA1  n  steps ize  at  wall 

N  Array  subscript 

The  following  variables  used  in  subroutine  PRFILE  are  described 
in  the  section  Description  of  Variables  in  Common:  DN,  EPSPL,  ETAINF, 

FC,  FCN,  FCP,  FI,  FIN,  FINN,  F2,  F2NN,  IE,  IM,  PNC,  VW,  XKETA,  XN,  Y, 
YOVDEL,  YOVTHT. 

Subroutine  READIN.  Subroutine  READIN  provides  the  input  of 
data  for  the  program.  Input  data  are  from  a  single  data  set  (SUB  DATA) 
when  calculating  a  noniterative  solution,  two  sets  for  the  initial 
iterated  solution  (SUB  DATA  and  HESSCP  OUTPUT),  and  the  two  plus  a  re¬ 
start  tape  for  the  subsequent  iterations.  READIN  also  resets  the  input 
flags  to  their  default  value  when  they  are  set  to  0  in  the  input  data 
and  performs  some  pressure  calculations  for  the  iterated  solution. 

USAGE: 

CALL  READIN 

Description  of  Variables. 

.1  Array  subscript 

LST  Alphanumeric  variable,  assigned  value  of  'LAST'  in  a  DATA 

statement 


49 


All  other  variables  are  described  in  the  section  Description  of 
Variables  in  Common,  except  for  LABEL  and  LSTC  which  are  described  in 
the  section. 

Description  of  Input  Data. 

Subroutine  RESTAR.  Subroutine  RESTAR  restarts  program 
execution  at  a  predetermined  body  coordinate  location.  On  the  first 
iteration,  when  the  marching  integration  reaches  the  restart  location 
(IRSTRT),  a  tape  is  written  containing  all  the  data  needed  to  begin 
execution  on  the  next  global  iteration.  RESTAR  first  reads  in  the  data 
written  from  the  previous  global  iteration,  and  then  reads  in  the  pre¬ 
vious  displacement  body  coordinates.  Finally  the  restart  tape  is  re¬ 
written  using  the  data  calculated  by  the  current  execution  and  the  co¬ 
ordinates  of  the  new  displacement  body  are  written  out  at  the  end  of  the 
current  execution. 

Subroutine  SOLVE.  Subroutine  SOLVE  calculates  the  solu¬ 


tion  for  a  general  parabolic  partial  differential  equation  when  the 
equation  is  written  in  the  standard  form. 


+  ^l^n  +  ^2*  +  ^3  +  =  ® 


In  order  to  obtain  a  solution  for  the  above  equation,  values  of  <|>  and 
the  first  and  second  derivatives  of  <j>  with  respect  to  n  from  a  pre¬ 
vious  iteration  or  from  a  previous  x  station  are  required.  These  data 
are  passed  through  the  argument  list  as  the  arrays  W1 ,  WIN,  and  WINN. 
The  new  values  of  <p  and  the  first  and  second  derivatives  are  returned 
through  the  argument  list  as  the  arrays  W2,  W2N,  and  W2NN. 

USAGE: 

CALL  SOLVE  (WINN,  WIN,  W2NN ,  W2N,  W2,  El,  Fll,  CRNI) 
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Description  of  Variables. 

A  Matrix  Inversion  Coefficient 

B  Matrix  Inversion  Coefficient 

CC  Matrix  Inversion  Coefficient 

D  Matrix  Inversion  Coefficient 

E  Array  of  Matrix  Inversion  Coefficients 

El  Value  of  e  at  wall 

F  Array  of  Matrix  inversion  coefficients 

Fll  Value  of  F  at  wall 

KON  Array  subscript,  descending  order 

N  Array  subscript 

W1  Value  of  $  at  previous  station 

WIN  3<j>/3n  at  previous  station 

2  2 

WINN  3  <f>/3n  at  previous  station 

W2  New  value  of  <j> 

W2N  New  value  of  3<j>/3n 

W2NN  New  value  of  32<)>/3n^ 

The  following  variables  which  are  used  in  subroutine  SOLVE  are 
described  in  the  section  Description  of  Variables  in  Common:  Al,  A2,  A3, 
A4,  CRN I ,  DN,  DS,  IE,  IM,  XN. 

Function  TLU.  Function  TLU  searches  in  the  array  X  for  the  two 
values  which  bracket  XSTAR  and  linearly  interpolates  for  the  corresponding 
value  in  the  array  Z.  The  returned  value  is  TLU.  If  XSTAR  is  not  within 
the  range  of  the  array  X,  TLU  is  set  equal  to  zero  and  the  error  flag  is 
set  equal  to  one.  The  array  X  must  be  strictly  increasing. 
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USAGE: 

TLU  (NTABLE,  Z,  X,  XSTAR,  NFLAG)  where  NTABLE  is  the  dimension  of 
the  arrays  Z  and  X  and  NFLAG  is  the  error  flag. 

Subroutine  VELDAT.  Subroutine  VELDAT  nas  been  added  for 
the  express  purpose  of  writing  out  data  sets  for  the  Navier-Stokes  code. 
The  routine  writes  the  solution  of  the  velocity  profile  at  the  upstream 
boundary  of  the  Navier-Stokes  region  to  be  used  as  boundary  conditions 
for  the  APLNS  code. 

USAGE: 

CALL  VELDAT 

Subroutine  WRITE1.  Subroutine  WRITE1  writes  out  the  input 
data  and  the  initialized  data.  Integer  quantities  are  written  first  in 
the  same  order  as  they  appear  in  the  input.  The  next  data  written  are 
the  freestream  quantities.  The  transition  parameters  are  then  written, 
followed  by  a  group  of  variables,  most  of  which  are  from  the  program 
input  and  appear  in  approximately  the  same  order.  Next  are  written 
the  integer  arrays  which  control  the  values  of  x  at  which  a  solution  is 
obtained  and  stations  at  which  full  profile  data  are  to  be  written. 
Finally,  the  values  of  axial  distance,  radius,  surface  distance,  surface 
pressure,  and  nondimensional  edge  velocity  are  written. 

USAGE: 

CALL  WRITE 1 

Description  of  Variables. 

J  Array  subscript 

The  other  variables  which  appear  in  this  subroutine  are  described 

* 

in  the  sections  Description  of  Output  and  Description  of  Variables  in 
Common. 
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Subroutine  WRITE2.  Subroutine  WRITE2  provides  the  output 
of  computed  results  after  a  converged  solution  is  obtained.  The  output 
is  in  two  sections.  The  first  section  is  written  at  every  station  at 
which  a  solution  is  obtained.  The  second  section  is  written  for 
specified  values  of  x.  The  first  line  of  output  contains  the  values  of 
geometry  related  variables  (such  as  x,  z,  and  s).  the  pressure  gradient 
term  B.  the  number  of  iterations,  and  the  station  counter,  K.  The 
second  line  contains  the  nondimens ionalized  edge  conditions  (Ue ,  Pg, 
etc.),  the  skin  friction  in  transformed  coordinates,  the  matched  eddy 
viscosity,  and  the  cumulative  total  of  the  number  of  iterations.  The 
third  line  gives  the  first  and  second  derivatives  of  the  velocity  at  the 
wall  and  at  the  second  grid  point  from  the  wall.  The  remainder  of  the 
first  section  is  omitted  at  x=0.  The  next  line  gives  transition  param¬ 
eters  and  the  following  line  contains  the  Reynolds  number.  In  the 
remaining  lines  in  the  first  section  are  skin  friction  coefficients,  the 
boundary-layer  thicknesses,  and  quotients  of  boundary-layer  parameters. 
In  the  second  section  of  the  output,  some  of  the  variables  which  are 
functions  of  the  normal  coordinate  are  written.  The  output  of  this 
section  contains  n.  y,  y/e,  y/6,  F(=U/Ue) ,  3F/9n»  y+,  u+,  Udef,  and  e+. 
USAGE: 

CALL  WRITE2 

Description  of  Variables. 

II  Counter  or  subscript  for  array  I  PR 

IPRNT  Counter  or  subscript  for  array  IPRFL 

The  values  of  n  and  the  functions  of  n  are  written 


KL 
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at  every  KLth  grid  point 
XOLD  Previous  value  of  x 

The  above  variables  are  in  common  and  are  also  described  in  the 
section  Description  of  Variables  in  Common.  The  other  variables  which 
appear  in  this  subroutine  (except  for  N,  an  array  subscript)  are 
described  in  the  section  Description  of  Output. 

Subroutine  ZRO.  Subroutine  ZRO  uses  subroutine  INTERS  to 
interpolate  in  the  arrays  ZA  and  RZ  for  the  values  Z  and  RO  corresponding 
to  the  value  X2  in  the  array  XSTA.  The  arrays  RZ,  XSTA,  and  ZA  contain 
values  of  radius,  surface  distance,  and  axial  distance  respectively  for 
an  axi symmetric  geometry. 

USAGE: 

CALL  ZRO  (JJ,  X2,  Z,  RO) 

Description  of  Input  Data. 

There  are  two  sets  of  input  files  to  be  used  in  obtaining  the 
iterated  solution,  i.e.,  SUB  DATA,  and  HESSCP  OUTPUT.  SUB  DATA  contains 
the  basic  data  for  execution,  and  HESSCP  OUTPUT  contains  the  Cp  distri¬ 
bution  at  the  edge  of  the  displaced  body.  SUB  DATA  will  be  described 
first. 

CARD  1  IREAD,  IRSRT,  ITER,  KHESS  (6(8X, 12) ) 

IREAD  Restart  flag,  0  no  restart  tape 

to  be  read  (or  written),  1  read 
(and  write)  restart  tape 

IRSRT  Station  number  for  restart  tape 

to  be  written 


ITER 


Iteration  counter 


CARDS  2-6 


CARD  7 


CARD  8 


CARD  9 


KHESS  Flag  for  using  Cp  distribution; 

0  do  not  read  in  Cp  data  from 
HESS,  1  read  in  Cp  data 
HEDR,  CASE,  ICARD2,  ICARD3,  ICARD4,  NN,  IBDN 
( 1 0A6 ,  2x,  A6/I3/I5/lx,I4,I5/9x,Il ) 

Standard  HESS  input  parameters,  described  in  HESS  PARMDATA 


and  HESSCOOR  DATA 
BETAS,  ROMAX 
BETAS 


ROMAX 

LABEL 


(6F12.6) 

Minimum  value  of  BETA,  the 
pressure  gradient  term.  If 
BETA  goes  lower  than  BETAS,  BETA 
is  set  equal  to  BETAS. 

Maximum  body  radius 
(18A4) 


Alphanumeric  data  for  identification  of  the  case.  Printed 
on  the  first  page  of  output. 


NOSE,  LAMTRB,  KVSAW,  KTRANS,  KTRNSN,  KONSET  ( 3 (8 x , 1 2 ) , 

3(7x,I3)) 

NOSE  Indicator  for  geometry;  1  for 

a  blunt  nose,  2  for  a  sharp 
nose  body.  A  value  of  0  is 
reset  to  1 

LAMTRB  Indicator;  1  for  laminar  flow, 

2  for  turbulent  flow.  A  value 


of  0  is  reset  to  1 
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KVSLAW 

KTRANS 

KTRNSN 

KONSET 

CARD  10  NIT1,  NIT2,  NIT3,  NC 

NIT1 

NIT2 

NIT3 

NC 

CARD  11  IE,  KEND,  UMAX 

IE 


Indicator:  0  for  Reichardt 
Eddy  Viscosity,  1  for  Van  Driest 
Eddy  Viscosity 

Indicator  for  transition  regime; 

0  for  r.o  transition  or  in¬ 
stantaneous  transition,  1  for 
transition  regime  (LAMTRB  must  be 
1) 

Subscript  for  XSTA  array  giving 
value  of  x  at  which  instantaneous 
transition  occurs 

Subscript  for  XSTA  array  giving 
value  of  x  at  which  transition 
regime  begins 

(4(8X,I2) ) 

After  solution  converges.  Ax 
is  doubled  if  NIT  is  less  than 
or  equal  to  0.  A  value  of  0  is 
reset  equal  to  3.  A  value  less 
than  zero  is  reset  equal  to  0. 
Recommended  value  is  5 

After  solution  converges.  Ax 
is  halved  if  NIT  is  greater 
than  or  equal  to  NIT2.  A  value 
of  0  is  reset  to  6.  A  recommended 
value  is  10 

If  NIT  =  NIT3,  Ax  is  halved,  X 
and  XI  are  recomputed  and  the 
iteration  loop  is  restarted  with 
the  smaller  values  of  X  and  XI. 

A  value  of  0  is  reset  to  9.  A 
recommended  value  is  20 

Specific  convergence  test;  0  for 
test  on  derivatives  at  wall,  1 
for  test  on  functions  at  all 
points  in  profile.  Recommended 
value  is  1 

(4(7X,I3) ) 

Number  of  grid  points  in  ETA  di¬ 
rection.  (program  dimensioned 
for  101,  the  recommended  value) 


■ . . -  -  -  — AA* 


_ 


_ _ _ 
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KEND 

UMAX 

CARD  12  KADETA,  KL,  IPFL,  IPRINT, 

IWT 

KADETA 

KL 

IPFL 

IWT 

CARD  6  CHI  CRT,  XBAR,  ATR 

CHI CRT 

XBAR 


Upper  index  on  K  (the  station 
counter).  Reconmended  value, 
300-600 

Number  of  values  in  geometry 
and  pressure  arrays.  Arrays 
dimensioned  for  500  values, 
recommend  no  more  than  101. 

If  more  than  101  points  are 
used,  changes  will  need  to  be 
made  to  the  dimensions  of  ar¬ 
rays  I  PR  and  IPRFL 


(5(8X, 12) ) 

Indicator;  if  zero,  ETAINF  is 
held  constant,  if  nonzero, 

ETAINF  can  be  increased  or  de¬ 
creased.  A  nonzero  value  is 
recommended 

Index  for  profile  output  DO 
loop;  1  gives  output  for  each 
grid  point,  2  for  every  second 
grid  point 

Number  of  values  in  IPRFL  array 
(maximum  100) 

Increment  in  station  number  at 
which  full  profiles  are  printed 

(3F12.6) 

Value  of  CHIMAX  at  which  transi¬ 
tion  regime  calculations  begin 
(arbitrary  if  KTRANS  =  0).  Ap¬ 
propriate  values  have  been  found 
to  be  between  2000  and  4000 

Measure  of  relative  length  of 
transition  region  (i.e.  Xend 

transition  "  XBAR*Xbeginning  of 
transition) 

For  Reichardt  Law,  XBAR  should 
be  about  65  percent  greater  than 
value  used  with  Van  Driest  Law. 
Recommended  values  are  approximately 
3.3  and  2.0,  respectively 
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ATR  Transition  constant  =  0.412 

CARD  14  UFS,  REINF,  PFS,  PSTAG,  (6F12.6) 

SCF,  RNOSE 


UFS 

REINF 

PFS 

PSTAG 

SCF 

RNOSE 

CARD  15  XK1,  XK2,  CONVRG,  ADTEST 

XK1 

X 

CONVRG 

ADTEST 


CARD  16  DX,  CRN  I ,  XKETZ,  ETAINF, 

REFLEN 

DX 

CRN  I 


Freestream  velocity,  ft/sec 

Freestream  Reynolds  number 
per  unit  length 

Freestream  static  pressure, 
psfa 

Freestream  stagnation  pressure, 
psfa 

Scale  factor  to  convert  unit 
length  for  REINF  (1.0  if  REINF/ft, 
12.0  if  REINF/in.  etc.) 

Nose  radius,  ft 

(6F12.6) 

Constant  for  Van  Driest  Inner 
Eddy  Viscosity  Law.  A  value 
of  0.0  is  reset  to  0.4,  the  com¬ 
mon  value 

Constant  for  Outer  Eddy  Viscosity 
Law.  A  value  of  0.0  is  reset  to 
0.0168,  the  common  value 

Convergence  criterion.  A  value 
of  0.0  is  reset  to  0.001,  the 
recommended  value 

Criterion  for  increasing  ETAINF 
(F(IE))-F( IE-4 )  is  compared  with 
ADTEST.)  A  value  of  0.0  is  re¬ 
set  to  0.001,  the  recommended 
value 

(6F12.6) 


Initial  x  stepsize 

Selector  for  finite  difference 
scheme;  1.0  for  fully  implicit 


♦ 
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XKETA 

ETAINF 

REFLEN 

CARD  17  XJFAC,  XKFAC,  TVC 

XJFAC 

XKFAC 

TVC 

CARD  18  DXMAX,  XLOCA,  DYDX, 

DX,  NOPRFS 

DXMAX 

XLOCAT 

DYDX 

DX 


(recommended),  0.5  for  Crank- 
Nicolson,  and  0.0  for  fully 
explicit  (not  recommended) 

Parameter  which  controls  the 
grid  spacing  in  the  ETA  direc¬ 
tion  (1.0  gives  an  equally  spaced 
grid,  a  value  greater  than  1.0 
gives  a  finer  grid  at  the  wall 
than  at  the  outer  edge.)  A 
value  of  0.0  is  reset  to  1.09, 
the  value  recommended  for  IE  =  101 

Maximum  value  of  ETA.  A  value  of 
0.0  is  reset  to  100  which  is  the 
recommended  value  for  a  flow 
which  is  initially  turbulent. 

For  an  initially  laminar  flow 
this  value  should  be  set  to  6.0) 

An  arbitrary  reference  length 

(2F12.6,4X,A3) 

Indicator;  0.0  for  2-D  flow, 

1.0  for  axially  symmetric  flow 

indicator;  0.0  for  body  only 
solution,  1.0  for  body  and 
wake  solution 

Indicator;  YES  to  include  TVC 
effects,  NO  without  TVC 


Maximum  x  stepsize  permitted  in 
the  calculations 

Location  on  the  body  of  the  up¬ 
stream  edge  of  the  Navier-Stokes 
region 

Slope  of  straight  taper  fit  of 
tail  of  body 

x  stepsize  for  use  in  the  Navier 
Stokes  code 


NOPRFS 


Number  of  profiles  on  the  body  in 
the  Navier-Stokes  region 


59 


CARDS  19a,  b,  c. . . 


(IPR(J) ,J  =  1,  IPRINT)  (1415) 


I  PR 


CARDS  20a,  b,  c... 


Array  of  subscript  for  array 
XSTA  giving  values  of  x  at  which 
complete  solutions  (with  pro¬ 
files)  are  printed  (array  di¬ 
mensioned  for  100  values) 

(IPRFL(J),  J  =  1,  IPFL)  (1415) 


IPRFL 


CARDS  21-21+IIMAX 

ZA 

XSTA 

RZ 

PZ 


LSTC 


Array  of  subscripts  for  array 
XSTA  setting  values  of  x  at 
which  complete  solutions  (with 
profiles)  are  printed  (array 
dimensioned  for  100  values). 

Also  at  these  stations  are 
printed  the  coordinates  of  the 
displaced  body.  Each  value  of 
IPRFL  must  also  be  a  value  of 
IPR 

ZA(J) ,  XSTA(J),  RZ (Z ) ,  PZ(J), 

LSTC  (4F12.6.A4) 

Array  of  coordinate  values  along 
the  body  axis 

Array  of  coordinate  values  along 
the  body  surface  (same  as  ZA  for 
a  flat  plate) 

Array  of  radius  values  for  a  body 

Array  of  pressure  ratio  values  - 
P/P0  where  PQ  is  the  local  stag¬ 
nation  pressure.  (These  are 
always  0.0  for  the  iterated  solu¬ 
tion  and  are  given  values  only 
when  KHESS  =  0) 

Alphanumeric  variable  used  to 
indicate  the  last  data  card. 

LAST  is  coded  on  the  final  card. 


See  Appendix  B,  Table  2.1  for  a  sample  SUB  DATA  data  set. 
Description  of  HESSCP  OUTPUT. 


CARD1-IIMAX 


X(J),  CP(J),  J  =  1 ,  UMAX 
(2F12.6) 
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X  Axial  coordinate  of  body  point, 

ft 

CP  Pressure  coefficient  at  body 

point 

See  Appendix  B,  Table  2.2  for  a  sample  HESSCP  OUTPUT  data  set. 

Description  of  Output  Data. 

ICBLINT  outputs  three  data  sets.  The  first  is  the  complete 
boundary- layer  solution  over  the  body  (BLAYER  OUTPUT).  The  second 
is  the  velocity  profile  at  the  upstream  edge  of  the  Navier-Stokes  region 
with  information  on  the  profile  included  (APLNS  VELDATA).  The  third  set 
is  the  coordinates  of  the  displaced  body  (B0DYC00R  DATA). 

This  section  contains  a  description  of  BLAYER  OUTPUT  in  three  sub¬ 
sections.  The  first  subsection  describes  the  first  section  of  the  out¬ 
put  which  is  the  input  and  initialized  data.  The  second  subsection 
describes  the  data  which  are  written  for  each  station  at  which  solutions 
are  obtained.  The  third  subsection  describes  the  miscellaneous  messages 
and  data  written  by  the  program.  The  output  from  a  sample  calculation 
of  turbulent  flow  over  a  cylinder  is  included  at  the  end  of  this  section 
(see  Appendix  C,  Table  1.1). 

Input  and  Initialized  Data. 

LABEL  ,  Alphanumeric  Data  for  identifi¬ 

cation  of  the  case  (written  by 
Subroutine  READIN) 

The  data  in  the  remainder  of  the  first  block  is  written  by 
subroutine  WRITE1  and  corresponds  to  the  input  data  with  the  following 


exceptions: 

Line  6 

RHOFS 

= 

P  ,  slugs/ft3 

00 

Line  10 

EPSVD 

= 

eVD  =  {uref/(p»U»L)]H’  where  L 

is  the  unit  length  in  REINF 

Next  in  the  program  output  are  the  arrays  IPR  and  IPRFL  (which  are 
described  in  the  section  "Description  of  Input").  Following  the  arrays 
IPR  and  IPRFL  are  the  arrays  ZA,  RZ,  XSTA,  PZ,  and  UEZ  with  the  array 
subscripts.  PZ  in  the  output  differs  from  the  values  of  PZ  in  the  input 
in  that  for  input  PZ  is  P/PQ  and  is  changed  by  the  program  to  P/Pref. 

UEZ  is  the  array  of  nondimensional  edge  velocities  calculated  by  the 
program  (Ug/U^).  For  a  flat  plate,  the  values  of  UEZ  are  unity. 

Station  Data. 

The  station  data  are  divided  into  two  groups.  The  first  group  of 
data  is  written  for  each  value  of  x  at  which  a  solution  is  obtained. 

The  second  group  —the  function  profiles— is  written  only  for  the  values 
of  x  corresponding  to  XSTA  ( IPFRL) .  In  the  first  group,  lines  5-13  are 


not  written 

if  x  = 

0.  In  the  sixth 

line,  REX  alone  is  written  if  KFS  =  0. 

If  KFS  f  0, 

REX  is 

written  followed 

by  the  wall  heat  transfer  (BTU/ft^/sec) 

and  two  heat 

transfer  coefficients; 

2 

one  in  BTU/in  -sec-°R,  the  other  in 

2 

lbm/in  -sec. 

GROUP  ONE 

Line  1 

S 

9 

Surface  Distance,  x 

XI 

= 

S,  Transformed  Surface  Distance 

Z 

9 

Axial  Distance,  z 

RO 

9 

Local  Radius,  r,  (Set  Equal  to 

One  for  a  Wedge  or  a  Flat  Plate) 

BETA 

= 

B  =  2c/Ue(dUe/dc) 

PP 

= 

dPE/ds 

NIT 

9 

Number  of  Iterations  at  Current 
Value  of  x 

K 

9 

Station  Counter 

DU  EDS 


dUE/ds 


Line  3 


Line  4 


Line  5 


Line  6 


CF 

9 

Skin  Friction  in  Transformed 
Coordinates  =  /e^ 

00 

EPSVD 

= 

eVD  =  {yref/  (p»u«>L)}J5‘  Where 

L  is  the  Unit  Length  in  REINF 

KEP 

9 

Subscript  in  EPSPL  Array  Where 
the  Inner  and  Outer  Eddy 
Viscosities  are  Matched  (Zero 
for  Laminar  Flow) 

NITTOT 

9 

Cumulative  Total  Number  of 
Iterations  in  x  Direction 

(F  =  U/Ue; 

Prime  Denotes  a/3n) 

F2N(1 ) 

= 

f; 

F2N (3 ) 

= 

F'  at  Two  Points  from  the  Wall 

F2NN ( 1 ) 

= 

F" 

w 

F2NN(3) 

= 

F"  at  Two  Points  from  the  Wall 

UPLUS 

= 

Friction  Velocity 

CHIMAX 

» 

Maximum  Value  of  Vorticity 
Reynolds  Number  -  x 

GAMMA 

» 

Transition  Intenmittency  Factor 
(Set  to  Zero  for  Laminar  Flow, 
Set  to  One  for  Turbulent  Flow) 

XI  BAR 

Transition  Intermittency 
Factor  Coordinate;  Initialized 
as  Zero,  Calculated  after 
Transition  Regime  is  Entered 

REX 

* 

Local  Reynolds  Number  Based  on 
Edge  Conditions 
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WALL 

HEAT 

TRANSFER  AND  HEAT  TRANSFER 

Line 

7 

CFE 

= 

CFINF 

= 

CF (TOTAL) 

= 

CFBAR*SQRT (REX) 

= 

CFINF*SQRT (REX) 

= 

Line 

8 

CF INF/2 

= 

Line 

9 

X/REFLEN 

= 

Z/REFLEN 

= 

RO/REFLEN 

= 

DELTA/X 

= 

DELSTR/X 

= 

Line 

10 

THETA/REFLEN 

= 

DELSTR/REFLEN 

= 

DELTA/REFLEN 

= 

DELSTRAXI/ 

REFLEN 

= 

Line 

11 

THETA/DELTA 

= 

DELSTR/DELTA 

= 

DELSTR/THETA 

= 

DELSTRAXI 

= 

Line 

12 

DELSTRK 

_ 

COEFFICIENTS  (If  KFS/O) 

Cf  =  2t  Jo  Ail 
f  wee 
e 

C,  =  2r  Jp  U2 

y  oo  oo 

00  /X 

Cf  =  X  J  Cf  dt 

Cf  /REX 
Cf  /REX 

00 

Cf  /2 

00 

x/reference  length 
z/reference  length 
r /reference  length 
6/x 
6*/x 

e/reference  length 
6*/reference  length 
6/reference  length 
6*  ./reference  length 

e/6 
6*/6 
6*/ 0 

6*  •  =  r  {/l+26*/r  -1 } 

aXl 

(Zero  for  2-D  Flow) 

6j*»  Incompressible  Boundary- 

Layer  Displacement  Thickness 
Computed  Only  for  Turbulent  Flow 
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DELSTR 

DELTA 

THETA 

RETHETA 

GROUP  TWO  PROFILES 

Block  1  ETA  , 

Y 

Y/THETA 

Y/DELTA 

F=U/UE 

FP(N) 

YPLUS 

U/UPLUS 

UDEF 

EPS+ 

PITOT 

N 

Miscellaneous  Messages  and  Data 


6*,  Boundary-Layer  Displacement 
Thickness 

6,  Boundary-Layer  Thickness 
(Value  of  y  where  U/Ug  *  0.995) 

e.  Boundary- Layer  Momentum 
Thickness 

Re.,  Edge  Unit  Reynolds 

0 

Number  Times  e 

n.  Transformed  Normal  Coordinate 

y.  Corresponding  Physical 
Coordinate 

y/e  (Zero  for  x  =  0) 
y/6  (Zero  for  x  =  0) 
Non-dimensional  Velocity 

II 

3n 

Law  of  the  Wall  Normal 
Coordinate,  yU+/v 

U/U+  Where  U+  =  (t^p)*5 

Defect  Velocity,  U"U°° 

U+ 

Eddy  Viscosity  (Zero  for 
Laminar  Flow) 

The  local  pitot  pressure  in 
p.s .t. 

Grid  Point  Number 


Subroutine  READIN  writes  the  message  "RESTART  TAPE  WRITTEN. 

K  =  #  KTPW  -  #"  after  the  station  data  has  been  written  on  the  restart 
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tape.  K  Is  the  corrent  value  of  K  and  KTPW  is  the  next  value  of  I  at 
which  the  restart  tape  will  be  written. 

For  transition  regime  calculations,  subroutine  EFFMU  writes  the 
message  "TRANSITION  BEGINS"  at  the  beginning  of  transition  and  the  message 
"TRANSITION  ENDS"  when  XI BAR  first  exceeds  2.0  with  the  Reichardt  Eddy 
Viscosity  Law  or  when  XIBAR  first  exceeds  4.0  with  the  Van  Driest  Eddy 
Viscosity  Law.  The  message  "TRANSITION  ENDS"  was  found  to  be  useful  when 
the  Van  Driest  Eddy  Viscosity  Law  was  used  in  calculations  but  less  useful 
with  the  Reichardt  Eddy  Viscosity  Law.  At  best,  the  message  indicates 
only  approximately  where  transition  ends. 

For  the  instantaneous  transition  model,  subroutine  CHANGE  writes 
the  message  "TRANSITION  BEGINS  INSTANTANEOUSLY"  and  the  values  of  X  and 
DX  when  transition  is  initiated. 

If  the  iteration  counter  NIT  reaches  the  value  of  NIT3,  the  stepsize 
is  halved,  new  values  of  X  and  XI  and  new  edge  properties  are  calculated 
and  the  program  attempts  to  obtain  a  converged  solution  for  the  smaller 
value  of  X.  In  this  case,  MAIN  writes  the  values  of  NIT,  DX,  X,  F2N(1), 
F2N1 ,  and  DIF. 

When  ETAINF  is  increased,  subroutine  ADDETA  writes  the  message 
"INTERMEDIATE  PROFILE  DATA— ETAINF  INCREASED  ISTOP  =  and  the  values 
of  DX,  XKETZ,  ETAINF,  and  ADTEST. 

If  ISTOP  exceeds  100,  the  message  "STOP  ***  ISTOP. GT. 100"  is 
written  by  MAIN  and  the  program  is  stopped.  If  after  a  converged 
solution  is  obtained,  the  derivative  of  F  at  the  wall  is  negative, 
subroutine  WRITE2  writes  the  message  "PROBLEM  TERMINATED.  NEGATIVE 
DF/DETA  INDICATES  THAT  THE  BOUNDARY  LAYER  HAS  SEPARATED"  and  stops  the 


program. 
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Normal  termination  of  the  program  is  indicated  by  the  message 
"THE  END  X  =  XSTA(IIMAX)"  or  the  message  "THE  END  K  =  KEND"  which  is 
written  twice  by  MAIN. 

In  the  above  descriptions,  all  quantities  which  have  the  apparent 
dimension  of  length  have  been  nondimensional ized  by  the  per  unit  length 
in  REINF. 

APLNS  VELDATA  contains  the  needed  parts  of  the  boundary-layer 
solution  at  the  beginning  of  the  Navier-Stokes  region.  All  information 
in  it  is  found  in  BLAYER  OUTPUT  and  is  written  again  for  convenience  in 
using  the  APLNS  program. 

BODY  COORDATA  contains  the  coordinates  of  the  displaced  body 
calculated  during  the  current  global  iteration.  This  data  set  is  used 
by  program  HESS  for  calculating  a  new  pressure  distribution  around  the 
body. 


See  Appendix  B,  Table  2.3  for  a  sample  BODY  COORDATA  data  set. 


DESCRIPTION  OF  PROCEDURE  DEVELOP 

After  the  viscous/ inviscid  solution  has  been  completed,  the  final 
displaced  body  coordinates  are  used  in  DEVELOP  EXEC  to  generate  the 
inviscid,  displaced  flowfield.  HESS  is  run  to  obtain  the  flow  quantities 
at  off  body  points  (u  &  v  components)  for  use  in  the  APLNS  program  input 
conditions.  Since  the  APLNS  code  requires  1800  points  (as  currently 
dimensioned)  and  HESS  is  capable  of  solving  only  400  points  per  execution, 
several  executions  of  HESS  are  required.  Each  execution  supplies  a  strip 
of  data  for  input  to  APLNS.  The  off-body  point  locations  are  developed 
and  the  data  set  written  out  in  DATFIELD  using  the  displaced  coordinates 
from  B0DYC00R  DATA.  DEVELOP  EXEC  controls  the  execution  of  FIELD  and 
HESSBLI  EXECs.  Execution  is  in  a  loop,  nominally  the  loop  is  executed 
five  (5)  times,  DATFIELD  is  executed  first,  then  with  the  data  set 
written,  HESSBLI  is  run,  writing  STREAM  DATA.  Subsequent  iterations  add 
data  to  STREAM  DATA  until  the  required  full  flowfield  is  completed 
(1800  points). 

See  Appendix  A,  Table  4.1  for  a  listing  of  DEVELOP. 

Program  DATFIELD 

Description  of  Program  DATFIELD.  FIELD  EXEC  controls  program 
DATFIELD.  DATFIELD  writes  out  the  data  sets  for  the  execution  of  HESS 
to  obtain  off-body  flow  conditions.  Input  data  include  the  final 
displaced  body  from  HESS/ICBLINT,  and  the  geometry  of  the  afterbody  from 
DATAMA.  DATFIELD  then  calculates  and  writes  the  off-body  points,  together 
with  the  displaced  body  coordinates  on  HESSCOOR  DATA,  the  HESS  data  set. 

See  Appendix  A,  Table  4.2  for  a  listing  of  DATFIELD  FORTRAN. 
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(See  Description  of  Output  for  a  list  of  variables  used  in  the  program.) 

Description  of  Input  and  Output. 

Input  Data.  DATFIELD  uses  two  input  sets:  INTER  DATA 
and  B0DYC00R  DATA.  INTER  DATA  has  been  described  in  Output  for  DATAMA. 
BODYCOOR  DATA  has  been  described  in  Output  of  ICBLINT. 

Output  Data.  DATFIELD' s  output  is  in  two  sets:  INTER 
DATA  and  HESS  COORDATA.  INTER  DATA  is  rewritten  to  show  the  iteration 
count  for  DEVELOP.  Only  the  counter,  IHESSI  is  changed.  HESS  COORDATA 
contains  the  HESS  data  set  with  the  off-body  points  where  a  solution  is 


desired. 

Card  1 

HEDR,  CASE 

(15A4,2x,A4) 

HEDR 

Alphanumeric  array  describing  the  body 
that  is  running 

CASE 

Alphanumeric  array  defining  the  case  run 

Card  2 

NB,NNU,NAXI,NCF, 

NOFF 

(511) 

NB 

Number  of  bodies  input 

NNU 

Number  of  nonuniform  flows 

NAXI 

Flag  for  axisymmetric  flows;  0  if  not 
axisymmetric,  1  if  axisymmetric 

NCF 

Flag  for  crossflow  case;  0  if  no  cross 
flows  included,  1  if  crossflows  included 

Card  3  blank 

NOFF 

Flag;  0  for  no  off-body  points,  1  for 
off-body  points 

Card  4 

IGEOMF,  IS IGF, 
ICURVN,  NONEWF, 
IFORMT,  NN1 

IGEOMF 


Flag  for  curved  elements;  0  for  no 
curvature,  1  for  curvature 
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Card  5 

Cards  6-NN 

Card  NN+7 


Card  NN+8 

Card  NN+8- 


ISIGF  Flag  for  source  densities;  set  to  1; 

pi ecewise-1 inear 

ICURVN  Flag  for  computation  of  element 

curvatures;  0  for  input  curvatures, 

1  for  internal  computation 

NONEWF  Number  of  new  crossflows  for  the  same 

body 

IFORMT  Flag;  if  0,  read  first  the  x,  then  the 

y  coordinates;  if  1,  read  x,y  coordinates 
in  2F12.6  format,  if  2  read  x,y  co¬ 
ordinates  in  F12.6,12x,F12.6  format 

NN1  Number  of  off-body  points 

IBDN 

IBDN  Body  sequence  number 

X,Y 

X,Y  Body  coordinates 


IGEOMF, 

rCURVN, 

IFORMT, 

ISIGF, 

NONEWF, 

NN1 

IGEOMF 

Flag  for  curved  elements;  0  for  no 
curvature;  1  for  curvature 

IS  IGF 

Described  above 

ICURVN 

Described  above 

NONL vF 

Described  above 

IFORMT 

Described  above 

NN1 

Described  above 

I  IBDN 

I  IBDN 

Body  sequence  number;  0  if  off-body 
points  follow 

XB,  YB 

XB,  YB  ,  Coordinates  of  off-body  points 
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See  Appendix  B,  Table  3.1  for  a  sample  HESS  COORDATA  data  set. 
Program  APLNS 

Description  of  Program  APLNS.  Program  APLNS  solves  the  complete 
Navier-Stokes  equations  using  an  Alternating  Difference  Implicit  (ADI) 
scheme.  Surrounding  boundary  conditions  are  required  to  initiate  execu¬ 
tion.  The  boundary  conditions  at  the  upstream  boundary  are  held  con¬ 
stant  while  the  remainder  of  the  flowfield  is  iterated  upon.  The  iter¬ 
ation  scheme  is  discussed  in  more  detail  later.  Due  to  the  large  grid 
size  used  within  the  program,  there  are  very  few  points  within  the 
boundary-layer  region.  Therefore,  the  results,  particularly  on  the  body, 
must  be  viewed  with  caution.  To  increase  the  accuracy  of  the  results 
obtained  from  the  APLNS  program,  the  user  is  advised  that  the  size  of 
the  region  of  interest  be  enlarged,  i.e.,  the  array  sizes  be  increased 
from  30X60  to  some  larger  values.  The  effect  of  this  will  be  to  make 
the  far  wake  boundary  conditions  used  in  the  program  more  appropriate  the 
larger  the  region  size,  the  better  these  conditions  become.  However, 
execution  time  will  increase  with  such  an  enlargement.  Care  must  be 
taken  with  regard  to  the  initial  flowfield  which  is  input  to  the  program. 
The  stepsizes  used  must  be  compatible  with  the  following  criterion: 

1.  There  must  be  at  least  three  or  four  points  within  the 
boundary- layer.  Therefore  AY  must  be  appropriatly  small. 

2.  The  ratio  aY/aX  must  match  the  slope  of  the  body  at  the  tail 
region.  Since  the  full  Navier-Stokes  equations  are  solved, 
the  solution  is  very  much  dependent  on  the  time  stepsize. 

For  high  velocities  a  very  small  time  differential 
(0.1  X  10  sec)  must  be  used. 


Nunerical  Procedure 


1.  Set  boundary  conditions  for  c^n+^  and  <|/n+^.  Coefficients 

for  u(n+Js)  and  v^n+J^  are  extrapolated  from  the  previous  time 
step.  Boundary  conditions  are  known  except  for  which 

must  be  Iterated  upon. 

2.  Solve  the  Navler-Stokes  equation  using  an  Alternating  Dif¬ 
ference  Implicit  (ADI)  scheme  to  get  s^n+^. 

3.  Solve  for  ij/n+^  using  a  direct  solver  method  with 
for  the  right  side. 

4.  Calculate  new  values  of  ij/n+^  and  ^n+^. 

5.  Iterate  on  and  £  until  au/ay)w  -  3u/3y  meets  the  convergence 
criterion. 

6.  Compute  u^n+1^  and  v^n+1^  from  the  new  stream  function. 

7.  Modify  the  stepsize.  If  necessary,  to  aid  convergence. 

Figures  10  and  11  show  the  flow  of  logic  and  the  calling  sequence 


of  subroutines  in  program  APLNS. 

See  Appendix  A,  Table  5.1  for  a  listing  of  program  APLNS. 
Description  of  Variables  in  Common. 


BLOCK 

VARIABLE 

DESCRIPTION  OF  USE 

REALS 

DFRODR 

1/p  3Fx/3r  lb-ft2/slug 

REALS 

DT 

At;  time  increment 
(sec) 

REALS 

DTA 

At;  time  increment 
(sec) 

REALS 

DT2 

At/ 2 

REALS 

DX 

Axial  length  increment 
(ft) 

A 
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REALS 

DY 

Radial  length  increment 
(ft) 

REALS 

DY2 

(ay)2  (ft)2 

ARRAYS 

EPA 

e,  epsilon  at  each 
axial  li cation  (eddy 
viscosity) 

REALS 

EPO 

Beginning  value  of 
epsilon  (ftVsec) 

REALS 

EPS 

Value  of  epsilon,  eddy 
viscosity,  if  <0 
e=0.01686Ue(ft2/sec) 

REALS 

FX 

Fx;  thrust  magnitude 
(lb) 

REALS 

FXRHO 

Fx/p  (lb-ft3/slug) 

INTEGS 

IDIMY 

Number  of  columns  for 
DMT010 

INTEGS 

10 

Input-output  unit  for 
any  saved  data 

INTEGS 

JCI 

Point  number  on  wall 
at  the  front  of  the 
regi on 

INTEGS 

JOEL 

Boundary-layer  edge 
point  at  the  front  of 
the  region 

INTEGS 

JDIMY 

Number  of  rows  for 
DMT010 

INTEGS 

L 

Number  of  points  per 
profile 

IRRAYS 

LAR 

Array  which  defines 
body  shape 

LOGICS 

LSWH1 

Control  parameter  set 
to  FALSE  value 

INTEGS 

LI 

Height  of  the  propeller 
zone  in  number  of 
points 
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INTEGS 

M 

Number  of  profiles 

INTEGS 

Ml 

Number  of  profiles  on 
the  body 

INTEGS 

NER 

Number  of  iterations 
on  W0RK6 

INTEGS 

NTS 

Number  of  iterations 
for  solution  convergence 

LOGICS 

OFF 

FALSE 

LOGICS 

ON 

TRUE 

LOGICS 

PRNT 

Controls  printing  for 
saved  output 

ARRAYS 

PSI 

ij>;  stream  function 
(ft-fysec) 

ARRAYS 

PSIB 

<|i;  used  in  W0RK6 
(ft3/sec) 

ARRAYS 

PSIO 

ij>;  on  top  of  region 
(ft3/sec) 

INTEGS 

PSKP 

Controls  frequency  of 
DU1CAL  and  EPSCA6 
execution 

REALS 

RHO 

p;  density  (slugs/ft  ) 

REALS 

R1 

Modifies  steps ize 

LOGICS 

SAVE 

Control  for  saving 
output 

PLOTZ 

SKEY 

Key  for  contour  plots 

PLOTZ 

SL 

Variables  in  PIXER 

PLOTZ 

SL1 

Variables  in  PIXER 

PLOTZ 

SU 

Variables  in  PIXER 

LOGICS 

SWHFX 

Control  for  propeller 

REALS 

TDX 

2 AX  (ft) 

REALS 

TDY 

2aY  (ft) 
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REALS 

TOL 

Tolerence  for  con¬ 
vergence  of  W0RK6 

REALS 

TRAMP 

Modification  factor  of 
1/p  (3Fx/3r)  if 
necessary 

ARRAYS 

U 

Axial  velocity  (ft/sec) 

ARRAYS 

UA 

Axial  velocity  in 

WORK 6  (ft/sec) 

REALS 

UO 

Boundary-layer 
velocities  at  the  front 
of  the  region  (ft/sec) 

ARRAYS 

U1 

Axial  velocity  at  upper 
boundary  of  the  region 
(ft/ sec) 

REALS 

U2 

Ug;  edge  velocity 
(ft/sec) 

ARRAYS 

UlO 

Initial  axial  velocity 
at  upper  boundary  of 
the  region  (ft/sec) 

ARRAYS 

V 

Radial  velocity 
(ft/sec) 

ARRAYS 

VA 

Radial  velocity  in 

W0RK6  (ft/sec) 

REALS 

X 

Axial  location  from 
front  of  the  region 
(ft) 

ARRAYS 

XI 

3u/3y  (sec-1) 

ARRAYS 

XIB 

3u/3y  calculation  at 
the  wall  from  stream 
function  in  W0RK6 
(sec-') 

ARRAYS 

XIW 

3u/3y  at  the  wall 
(sec"1 ) 

REALS 


Y 


Radial  distance  from 
body  centerline  (ft) 


Description  of  Subroutines. 


DRFAX6.  DRFAX6  is  the  main  driver  for  the  axi symmetric 
case.  The  routine  calls  the  major  routines  and  tests  for  convergence  of 
W0RK6  calculations.  The  routine  first  calls  two  subroutines  which  read 
in  data  and  set  up  initial  conditions.  Next,  the  routine  begins  the 
iteration  procedure  on  the  stream  function,  the  main  iteration  loop  of  the 
program.  After  the  convergence  is  met,  the  u  velocities  on  the  upper 
boundary  are  modified  every  PSKP  time  steps  and  modifications  to  the  eddy 
viscosity  are  made  every  PSKP+1  time  steps.  The  iteration  on  W0RK6  is 
repeated  NTS  time  steps  with  Dt  modified  depending  on  the  number  of 
iterations  required  for  convergence  of  the  previous  time  step.  If  more 
than  three  iterations  were  required,  Dt  is  decreased  by  R1 .  If  less  than 
three  iterations  are  required,  the  time  step  Dt  is  increased  by  the 
same  factor  R1 .  After  NTS  time  steps  or  the  number  of  iterations  in 
W0RK6  exceeds  NER,  DRFAX6  calls  a  routine  to  print  the  current  stream 
function,  au/ay  distribution,  and  the  u  and  v  velocity  components. 

Subroutine  DU1CAL.  DU1CAL  calculates  a  Au  that  is  added  to 
the  u  velocity  on  the  upper  boundary  of  the  region.  The  routine  is  called 
periodically  to  update  this  boundary  condition.  A  modified  value  for  the 
stream  function  is  used  to  calculate  the  Au.  The  routine  then  prints  out 
the  profile  number  (I),  the  modified  upper  boundary  condition  (UI),  the 
difference  between  the  current  velocity  and  the  initial  velocity  (DDU1), 
the  initial  velocity  (Ul 0) ,  and  the  initial  stream  function  value  (PSIO). 
It  is  important  to  note  that  only  velocities  on  the  upper  boundary  are 
altered  in  this  routine. 


1 


Li 
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Subroutine  EPSCA6.  EPSCA6  evaluates  the  eddy  viscosity 
law.  It  is  called  periodically  as  DU1CAL  is  called.  The  equation  used  is: 


i  =  C2b,  u  =  U/Ug,  r  =  r/b,  a2  =  2.0,  C2  =  0.035. 
The  routine  prints  out  the  following  values 


11 

Defined  above 

12 

Defined  above 

B 

Value  of  b  in  above  equation 

DB 

aB 

Dll 

aI1 

FI,  F2,  F3 

Terms  in  the  eddy  viscosity  law 

EPN 

New  value  of  e(x) 

EPA 

Original  value  of  e(x) 

These  updated  values  are  then  returned  to  the  program  and  used  to  modify 
the  flowfield  conditions. 

Subroutine  INTAL6.  INTAL6  is  responsible  for  calculating 
the  inviscid  stream  function  from  flowfield  velocities  read  in  from 
Unit  8  data.  A  simple  integration,  which  has  been  found  to  be  accurate. 
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is  used  to  obtain  the  inviscid  stream  function.  Since  the  velocities  at 
every  station  in  the  flowfield  are  known,  the  stream  function  is  calculated 
using  a  trapazoidal  rule: 


a*1j  ■  <ul/i  *  ui,j-ly1-2>ay/Z 

*1J  ‘  *1,0-1  +  S*1,0 

with  iji  =  0  at  solid  body  boundaries  and  for  zero  radius  in  the  wake. 

The  inviscid  flowfield  is  printed  out  giving  values  of  ip(PSI)  and  u  and 
v  velocity  components.  The  value  of  r  at  the  edge  of  the  boundary-layer 
at  the  front  of  the  region  ( Y(JDEL) )  is  printed  along  with  the  boundary- 
layer  thickness  (Y(ODEL)-Y(J(I))) ,  and  Ug  (U2).  The  boundary-layer 
velocity  profiles  are  read  in  from  Unit  11,  and  the  values  at  the  point 
in  question  Y(J)  are  interpolated  using  a  3  point  Lagrangian  formula. 

The  momentum  increase  (FX(1 ,L1 ) )  is  printed  out  giving  the  initial  values 

of  thrust  computed  from  the  following  equation; 

/6 


Fx/p  =  Fxn  =  2- 


SU 


0 


y(Ll)  DX 


5 —  ;  where  SUn 

£rw  U 


'o 


U(Ue  -  U)dy 


In  attempting  to  correct  for  tip  losses  the  following  expression  was  used: 
Fx, 


Fx/p 


'0  J 


,  .  ...  T  r-(Ll-5)Ar~ 
1  +  cos  4Ar 


r  for  (Ll-5)Ar<r<(Ll-l )Ar 


where  LI (dr)  is  the  height  of  the  propeller  region.  The  value  of  6 

(DELSR2)  is  calculated  by; 

✓  6 


*2 


*2 


(1  -  j}  )rdy 
e 


★ 

and  the  corresponding  value  of  e(EPS)  is  given  by:  e  =  0.01686  Ufi. 

The  boundary-layer  velocity  profile  (U0)  is  printed  out  along  with  the 
values  of  X  and  Y  for  the  region  in  question.  The  X  distances  are  taken 
from  the  front  of  the  region  and  the  Y  values  are  measured  from  the  axis 
of  the  body.  Finally  (1/p)  3Fx/sr(FXRH0)  is  printed  completing  the  printed 
output  from  the  routine. 

Subroutine  PUTIN6.  PUTIN6  reads  ir.  values  for  the  flags 
and  various  other  parameters  in  the  program  and  prints  them  out  in  namelist 
form.  The  routine  defines  the  body  shape  by  the  LAR  array  which  is  also 
printed  out.  The  initial  values  for  variables  which  control  the  key 
values  for  the  contour  plot  of  3u/3y(XI)  are  printed.  For  a  description 
of  these  parameters  see  the  section  entitled  Description  of  Input. 

Subroutine  W0RK6.  W0RK6  calculates  the  updated  values  of 
£  which  are  compared  in  DRFAX6  and  then  modified  if  necessary.  Since  this 
routine  essentially  solves  n  equations  in  n  unknowns,  only  the  numerical 
procedure  is  important  to  the  understanding  of  the  workings  of  W0RK6. 

The  use  of  W0RK6  is  best  seen  by  referring  to  the  section  entitled 
Numerical  Procedure,  steps  1-5. 

Description  of  Input. 


Input  Variables  on  Unit  5. 

Namelist 

Variable 

Description  of  Use 

X* 

HIST 

L 

Number  of  points  per 
profile,  suggested  value 
<30 

ILIST 

LI 

Number  of  points  on 
propeller,  suggested 
value  <12 

ILIST 

M 

Number  of  profiles, 
suggested  value  <60 
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ILIST 

Ml 

» 

Number  of  profiles  on 
body,  suggested  value 

U5 

HIST 

PSKP 

9 

Time  parameter  controlling 
frequency  of  DU1CAL  and 
EPSCA6,  suggested  value 

4-7 

ILIST 

NER 

9 

Number  of  iterations  on 
W0RK6,  suggested  value  40 

ILIST 

NTS 

9 

Number  of  time  steps  for 
convergence,  suggested 
value  35 

ILIST 

10 

9 

Unit  number  for  saved 
output,  suggested  value  9 

RLIST 

DT 

9 

Time  increment  (sec), 
suggested  value  1  x  10  - 

1  x  10-5 

RLIST 

RHO 

9 

Density  of  fluid  (slugs/ft^ 
suggested  value  7.535  x  10" 

RLIST 

EPS 

9 

Eddy  viscosity  (if  EPS<0, 
EPS  =  0.01 686*Ue) , 

suggested  value  0.0 

RLIST 

TOL 

9 

Tolerance  value  for 
convergence  of  W0RK6, 
suggested  value  0.02-0.05 

RLIST 

TRAMP 

9 

Modifying  factor  for 
propeller,  suggested 
value  0.1 

RLIST 

EPO 

9 

Beginning  value  of  epsilon, 
suggested  value  0.0 

RLIST 

SCALE1 

9 

Length  scaling  factor  for 
nondimens ional  HESS/ICBLINT 
1/0 

RLIST 

SCALE2 

9 

Velocity  scaling  factor 
for  nondimens ional 
HESS/ICBLINT  I/O 

LUST 


SWHFX 


T  No  actuator  disk 
F  Include  actuator  disk 
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(—  ,  -  .2. 


LUST 


SAVE 


T  Save  data  on  previous 
unit  defined  by  variable 
10 

F  Do  not  save  data 


See  Appendix  A,  Table  3.1  for  sample  APLNS  data. 
Input  Data  on  Unit  11 . 


Variable 

DYDX 

DY 

NUM 

DELTA 

UE 

YD 

UUE 


Description  of  Variable 
dy/dx;  slope  of  the  body  at  the  tail 
Ay;  nondimensional  radial  lingth  increment 
Number  of  boundary-layer  velocities  given 
Nondimensional  boundary-layer  thickness 
Nondimensional  edge  velocity 
y  locations  in  the  boundary  layer 
U/Ug;  normalized  boundary- layer  velocities 


See  Appendix  B,  Table  3.2  for  sample  APLNS  Unit  11  input  data. 
Input  Data  on  Unit  8. 

X  Nondimensional  axial  distance  from  the  front  of 

the  body 

Y  Nondimensional  radial  coordinate  from  the  body 

centerline 

U  Nondimensional  axial  velocity  component  (INVISCIt) 

V  Nondimensional  radial  velocity  component 

(INVISCID) 


See  Appendix  B,  Table  3.3  for  sample  APLNS  Unit  8  data. 
Description  of  Output. 

Output  Variables.  Listed  in  order  of  appearance. 


Origin 

Variable 

Description  of  Use 

PUTIN6 

L 

Number  of  points  per 
profile 
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PUTIN6 

LI 

Number  of  points  in  the 
propeller  zone  (height) 

PUTIN6 

M 

Number  of  profiles 

PUTIN6 

Ml 

Number  of  profiles  on  the 
body 

PUTIN6 

PSKP 

Control  variable  concerning 
frequency  of  DU1CAL  and 
EPSCA6 

PUTIN6 

NER 

Number  of  iterations  on 
W0RK6  (maximum) 

PUTIN6 

NTS 

Total  number  of  time  steps 
for  solution  convergence 

PUTIN6 

10 

Unit  number  for  saved 
output 

PUTIN6 

DT 

Time  increment  (s) 

PUTIN6 

RHO 

Density  (slugs/ft3) 

PUTIN6 

EPS 

Eddy  viscosity  (ft  /s) 

PUTIN6 

TOL 

Tolerance  value  for  con¬ 
vergence  of  W0RK6 

PUTIN6 

TRAMP 

Modification  factor  for 
derivative  of  thrust  in 
the  radial  direction 

PUTIN6 

EPO 

Eddy  viscosity  at  beginning 
of  problem  (ftz/s) 

PUTIN6 

SCALE 1 

Length  scaling  factor  (ft) 

PUTIN6 

SCALE2 

Velocity  scaling  factor 
(ft/s) 

PUTIN6 

SWHFX 

Control  parameter  for 
propeller 

PUTIN6 

SAVE 

Control  parameter  for 
printing  saved  output 

PUTIN6 

LAR 

Array  defining  body  points 
in  terms  of  grid  location 
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PUTIN6 

DX 

Axial  length  increment(ft) 

PUTIN6 

DY 

Radial  length  increment  ( ft) 

PUTIN6 

JOEL 

Point  number  at  the  edge  of 
the  boundary- layer 

INTAL6* 

PS  I 

Inviscid  stream  function 
(ft3/s) 

INTAL6* 

U 

Inviscid  axial  velocity 
(ft/s) 

INTAL6* 

V 

Inviscid  radial  velocity 
(ft/s) 

INTAL6 

Y(JDEL) 

Radial  distance  at  the 
edge  of  the  boundary- 
layer  (ft) 

INTAL6 

Y(JDEL)-Y(J(I)) 

Boundary-layer  thickness 
(ft) 

INT.AL6 

U2 

Axial  velocity  at  the  edge 
of  the  boundary-layer  (ft/s) 

INTAL6 

FX 

Force  distribution  of  the 
actuator  zone  (lb^) 

INTAL6 

DELSR2 

Square  of  the  boudary- 
layer  displacement 
thickness  (ft’) 

INTAL6 

EPS 

Eddy  viscosity  (ft2/s) 

INTAL6 

UO 

Boundary-layer  profile  (ft/s) 

INTAL6 

Y 

Radial  distances  of  grid 
points  from  body  centerline 
(ft) 

INTAL6 

X 

Axial  distances  of  grid 
points  from  front  of  region 
(ft) 

INTAL6 

FXRHO 

Derivative  of  propeller 
force  with  respect  to 

radial  direction 
(lbf-ft3/s1ug) 
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DRFAX6 

TIME 

Time  from  the  beginning 
of  the  procedure  (s) 

DRFAX6 

XI 

Derivative  of  axial 
velocity  with  respect  to 
radial  direction  (1/s) 

W0RK6 

N 

Dimension  of  working 
matrix  in  W0RK6 

WORKS 

IDIMY 

Dimension  of  working 
matrix  in  W0RK6 

W0RK6 

JDIMY 

Dimension  of  working 
matrix  in  W0RK6 

WORK6 

I 

Profile  number 

W0RK6 

J 

Point  number  on  profile 

WORK6 

XI(I,J) 

XI  at  point  I ,J  (1/s) 

WORK6 

XI W 

Vorticity  at  the  wall  of 
the  body  at  point  I,J  0/s) 

DRFAX6 

ITER 

Iteration  number  in  W0RK6 

DRFAX6 

XMAX 

Maximum  difference  between 
|(XI-X1W)/XI|  for  all  wall, 
points 

DRFAX6 

TS 

Time  step  number 

DU 1  CAL 

ui(i) 

Updated  axial  velocity  at 
the  top  of  the  region  (ft/s) 

DU 1  CAL 

DDU1 ( I ) 

Difference  between  updated 
and  original  upper  boundary 
axial  velocity  (ft/s) 

DU1CAL 

UIO(I) 

Initial  axial  velocity  at 
the  upper  boundary  (ft/s) 

DU1CAL 

PSIO 

Stream  function  at  the  top 
of  the  boundary  (ft3/s) 

EPSCA6 

B 

Mixing  zone  width  (ft) 

r  H 

EPSCA6 

11 

I,(x)  =  J  Qrdf 
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EPSCA6 

12 

J2(x>  -  J0  If  p  d? 

EPSCA6 

DB 

Change  in  the  value  of  the 
mixing  width  (ft) 

EPSCA6 

Dll 

Change  in  the  value  of 
I-|(x)  integral 

EPSCA6 

FI 

Term  in  the  eddy  viscosity 
law 

EPSCA6 

F2 

Term  in  the  eddy  viscosity 
law 

EPSCA6 

F3 

Term  in  the  eddy  viscosity 
law 

EPSCA6 

EPN 

New  value  of  the  eddy 
viscosity  at  profile  I 
(fWs) 

EPSCA6 

EPA 

Initial  value  of  the  eddy 
viscosity  for  profile  I 
(ft2/s) 

DRFAX6* 

PSI 

Final  stream  function 

(ft3/s) 

DRFAX6* 

XI 

Final  distribution  for 

XI  (1/s) 

DRFAX6 

V 

Final  values  of  radial 
velocities  (ft/s) 

DRFAX6 

U 

Final  values  of  axial 
velocities  (ft/s) 

*  _ 

Denotes  p  nted  within  mentioned 
routine  DKfOlO 

subroutine  through  a  call  to 

See  Appendix  B,  Table  3.4  for  sample  APLNS  OUTPUT  data. 
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a  TABLE  3.3  has  been  appended  from  180  to  85  lines. 
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